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ABSTRACT 

Earth-Crossing Asteroids (ECAs) are those asteroids whose orbit cross-section 
can intersect the capture cross section of the Earth as a result of secular gravitational 
perturbations. This thesis provides a framework for understanding the origin, nature, and 
types of ECAs. The change in velocity requirements to achieve a two Earth radii 
deflection for long and short-term warning scenarios are developed. Next, a method of 
developing hypothetical Earth colliding asteroid orbits is presented. These hypothetical 
orbits are used in two ways: (1) to evaluate the ability of Dance of the Planets, a solar 
system simulation model developed by Applied Research and Consulting, Inc., to 
accurately propagate orbits of imported asteroid orbits, and (2) to analyze the sensitivity 
of deflection distance to variation in deflection angle and orbital parameters of a given 
orbit. Inaccuracies during importation of data precluded the use of Dance of the Planets 
for the purpose of sensitivity analysis. The program does provide an excellent tool for 
visualization of ECA scenarios. Consequently, a simpler orbital model was developed to 
provide a Earth miss distance sensitivity analysis. With one asteroid orbital period 
warning the minimum change in velocity to deflect an asteroid two Earth radii is 
approximately 0.135 m/s and the optimal deflection is along the flight path. Maximum 
deflection occurs when the deflection is applied at perihelion. The miss distance decreases 
markedly with increase in true anomaly until it is a minimum at aphelion. 
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[. INTRODUCTION 


Impacts by Earth-crossing Asteroids (ECAs) and comets, collectively known as 
Near-Earth Objects (NEOs), pose a significant and unique challenge to the scientific 
community. The study of these space-borne bodies has spanned a broad range of 
disciplines: mechanical, electrical, aeronautical, and astronautical engineers, 
astronomers, nuclear physicists, to name but a few. Various camps have, at different 
times, proposed new search strategies and detection methods to ensure a proper tally of 
potential colliders, forwarded techniques, both nuclear and non-nuclear, to mitigate the 
disaster a colliding body may impart, and designed missions to intercept, rendezvous, and 
study a few of the closer bodies. The recent impact of comet Shoemaker-Levy 9 with 
Jupiter only served to stimulate additional interest in the subject of NEOs. 

Recognizing the potential seriousness of such events, the United States Congress 
in 1992 mandated that the National Aeronautics and Space Administration (NASA) 
conduct two workshops to study two major research areas of NEOs: Detection and 
Mitigation. The United States House of Representatives, in NASA Multiyear 
Authorization Act of 1990 said, in part: 

The committee believes that it is imperative that the detection rate of 
Earth-orbit-crossing asteroids must be increased substantially, and that the 
means to destroy or alter the orbits of asteroids when they threaten 
collision should be defined and agreed upon internationally. 

The chances of the Earth being struck by a large asteroid are extremely 
small, but since the consequences of such a collision are extremely large, 
the Committee believes it is only prudent to assess the nature of the threat 
and prepare to deal with it. We have the technology to detect such 
asteroids and to prevent their collision. 

The Committee therefore directs that NASA undertake two workshop 
studies. The first would define a program for dramatically increasing the 
detection rate of Earth-orbit-crossing asteroids; this study would address 
the costs, schedule technology, and equipment required for precise 
definition of the orbits of such bodies. The second study would define 
systems and technologies to alter the orbits of such asteroids or to destroy 
them if they should pose a danger to life on Earth. The Committee 



recommends international participation in these studies and suggest that 

they be conducted within a year of the passage of this legislation. [Ref. 1 ] 

As a result two conferences were conducted. The first, the NASA International 
Near-Earth-Object Detection Workshop, was conducted in three sessions from June 
through November 1991. Their work concentrated on improving the rate at which ECAs 
are discovered; the results are documented in reference 1. The second conference, The 
Near-Earth-Object Interception Workshop, was held in January 1992 and hosted by Los 
Alamos National Laboratory. It concentrated on the issues surround the mitigation of 
ECAs and associated technologies. Even more recently, an ECA conference was held in 
March 1995 at Lawrence Livermore National Laboratories. Conference proceedings have 
not yet been published. 

These conference reports address the majority of the issues surrounding NEOs; 
however, there is much work remaining. Advances in detection technology such as 
improved sensors, computer search algorithms, and space radar support would 
significantly speed the rate at which NEOs are identified. Space missions to candidate 
asteroids can provide essential information regarding the structure and composition of 
asteroid bodies. Optimization of intercept trajectories and continued improvements to the 
space launch vehicles will all be directly applicable to the mitigation of these potentially 
dangerous objects. 

At the onset of this project there were three main goals: (1) Conduct a thorough 
and exhaustive survey of the current literature in the area of ECAs, (2) Analyze and 
select from commercially available software applications the candidate most likely to aid 
in the visualization of asteroid mitigation through deflection, and (3) Develop mitigation 
scenarios and determine the sensitivity of the circular solution to elliptical orbits. 

Research under the Space Warfare Research Program, sponsored by the Air 
Force’s Space Warfare Center, was divided into four sections. Two major components of 
NEOs, ECA’s and Earth-crossing Comets (ECCs), are described and categorized. 
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Deflection options are presented and the maximum size asteroid which could be deflected 
is estimated. The most effective deflection angle is established and the required change 
in velocity to perturb an object by two Earth radii is determined. Methodology for 
developing hypothetical asteroid orbits is also presented. A commercially available 
software application. Dance of the Planets, is evaluated for its usefulness in verifying the 
deflection of assailant objects and testing the mathematical solutions for optimal 
deflection. Finally, an analysis of sensitivity of miss distance to variations in orbit shape 
and deflection direction is conducted. The issues compiled in this thesis address only a 
small portion of the asteroid mitigation problem, and an even smaller portion of the 
broader subject of ECAs. 
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II. NEAR EARTH OBJECTS 


A. INTRODUCTION 

There are two broad categories of objects whose orbits bring them close to that of 
the Earth: asteroids and comets. Astronomers classify the objects into one of the two 
categories based upon their telescopic appearance. Any object which appears to be star- 
like is called an asteroid. If the object has a visible atmosphere or a tail then it is a comet. 
[Ref. 1 ] This difference is, in part, due to the composition of the of the object. Asteroids 
have no atmosphere and may have physical and compositional properties ranging from 
loosely aggregated cometary ices to solid metal. [Ref. 2] Comet nuclei are composed of 
a complex mixture of volatile ices, water, and hydrocarbon and silicate grains. As a 
comet is heated by the Sun as it approaches perihelion, there is a noticeable out-gassing of 
evaporative material. [Ref. 3] Further discussions on each broad classification of Near- 
Earth Objects (NEOs) are provided below. 

B. ASTEROIDS 

When viewing the solar system, shown in Figure 1 Figure 1, the empty gap 
between Mars and Jupiter is readily apparent. Kepler contemplated a “missing planet” as 
he studied the solar system, as described by the measurements of Tycho Brahe. [Ref. 4] 



Figure 1. Depiction of Solar System Through Jupiter. From Ref. 4. 


Bode’s law, named after Johann Bode’s effort to stir interest in an unknown planet 
in 1772, uses a simple relationship to generate the mean distances in Astronomical Units 
(AU) of the principal planets. The relationship is generated by writing down the series 0, 
3, 6,12,..., add 4 to each number and dividing by 10. As shown in Table 1, Bode’s law 
generates fairly accurate locations for all the planets except Neptune and Pluto. 


PLANET 

BODE’S LAW 
DISTANCE 

ACTUAL 

DISTANCE 

Mercury 

0.4 

0.39 

Venus 

0.7 

0.72 

Earth 

1.0 

1.00 

Mars 

1.6 

1.52 

Asteroids (average) 

2.8 

2.65 

Jupiter 

5.2 

5.20 

Saturn 

10.0 

9.52 

Uranus 

19.6 

19.28 

Neptune 

38.8 

30.17 

Pluto 

77.2 

39.76 


Table 1. Bode’s Law vs. Mean Planetary Distance 




This fits the actual positions of the planets through Saturn very well, except for the 
position at 2.8 AU. This position, between the orbit of Mars and Jupiter remains empty, 
except for the asteroid belt. [Ref. 4] 

1. Early History 

In view of Bode’s law, a group of 24 astronomers formed a society in Europe in 
1800 to solve the problem of a planet, surmised to be missing, between Mars and Jupiter. 
Each astronomer was given a region of the zodiac. [Ref. 4], 

Giuseppe Piazzi, in Palermo, Italy, was already engaged in a star charting project. 
He located a dim, uncharted ‘star’ that shifted position from night to night. This was 
Ceres, discovered on January 01,1801. He observed it for 41 days before it was lost due 
to bad weather and the illness of Piazzi. [Ref. 4], 

Karl Friedrich Gauss became involved with the intricate mathematical problem of 
calculating an adequate recovery orbit for Ceres. With his help, Ceres was sighted on 
December 07, 1801. Gauss’ ingenuity was of great significance to the field of celestial 
mechanics. The basic elements of Ceres were eccentricity of 0.08, inclination of 11°, and 
semi-major axis of 2.77 AU. [Ref. 4]. 

Heinrich Wilhelm Olbers found a second unknown ‘planet’, Pallas, in March 
1802, while helping Gauss observe Ceres. Gauss next calculated an orbit for Pallas. 

Many scientists contributed in the calculation of the orbital elements, and in the 
development of perturbation theories to account for the influence of Jupiter. This was 
one of the first documented efforts at the development of planetary theory from which 
highly accurate planet ephemerides are possible. [Ref. 4], 

Discovery of additional asteroids progressed rapidly over the course of the next 
century, particularly with the implementation of photography in the search process in 
1891. By 1900,463 asteroids had been discovered, and by 1950, 1568. By 1993, 
approximately 5500 numbered asteroids had been documented. Approximately 300 



numbered asteroids are added annually. Numbered asteroids have orbits confirmed by 
observations at two or more oppositions (when the asteroid-Earth-Sun angle equals 180°). 
Many others asteroids have only provisional designations. [Ref. 4], 

2. Origin and Nature of Asteroids 

There are many plausible explanations for the existence and formation of 
asteroids. One scenario suggests the asteroids formed from planetesimals which were 
never able to accrete into planet sized bodies due to the disturbing influences of a proto- 
Jupiter. There are many inconsistencies with this theory. First, from the planetesimal 
model, there should be between one and two Earth masses of planetesimals in the region 
of the asteroid belt. All of the asteroids together, however, are no more than 0.08% of the 
mass of the Earth. Some mechanism is responsible for this loss of mass. Secondly, 
individual asteroids have high relative velocities. They do not fit the orderly picture of 
planetesimals analogous to Saturn’s rings: bodies lying in nearly the same plane with low 
relative velocities between neighboring bodies. Again, some mechanism is responsible 
for the high relative velocities which produce collision fragmentation and not collision- 
accretion. Finally, many stony and metallic meteorites, derived from asteroids, have clear 
indications of melting, elemental differentiation, and slow cooling, similar to that found 
in the planets. The process of this development is not fully understood. [Ref. 4], 

There has been much debate as to the origin of near-Earth asteroids. Because 
these are planet crossing asteroids, they have a dynamic lifetime (the average time before 
a planet crossing asteroid impacts one of the inner planets) of approximately 3 x 10 7 yr. 
As such, there must be some mechanism responsible for their replenishment. There are 
two main theories, divided largely between observers and dynamicists, which attempt to 
explain the source of these objects. Observers generally believe that the near-Earth 
asteroids were derived from the main asteroid belt due to spectropically similar objects in 
the main belt. Dynamicists believe they may be largely of cometary origin. They do not 



believe that there is sufficient dynamical interaction in the main asteroid belt to resupply 
the near-Earth orbit asteroid population. [Ref. 3] 

Asteroids are numbered in the order they are discovered with orbits verified with a 
minimum of three observations. The majority also have names assigned. The first 
asteroid discovered, 1 Ceres, is approximately 932 km in diameter, and represents about 
30% of the total mass of asteroids. 4 Vesta and 2 Pallas, approximately 528 km and 523 
km, respectively, are next in size. There are approximately 30 other main belt asteroids 
with diameters larger than 200 km. For objects with the same density, the mass varies as 
the third power of their linear dimensions. From this, it is easy to ascertain that most of 
the total mass of the asteroids is found in the few largest. Ceres is just below the limit for 
the class-3 satellite size category of 1000-1600 km, of which Saturn and Uranus each 
have four. [Ref. 4], 

3. Asteroid Close Approaches 

Although there are a large number of asteroids, their number is limited by our 
ability to detect them. This remains one of the biggest challenges in the field. Most 
newly discovered asteroids are dimmer and usually smaller than previously discovered 
asteroids. By far the majority are fragments of ancient parent bodies, unlike Ceres which 
has retained most of its original shape and mass. [Ref. 4], 

Near Earth Asteroids (NEA) are of interest because they lie well inside the main 
asteroid belt, and they do not have long term stability because they will eventually either 
1) collide with one of the inner planets, or 2) be ejected from this region by a near¬ 
collision encounter. Earth impact rate for objects of potentially cataclysmic size is 
estimated at 3-4 per million years. No such collisions have been yet been predicted, 
however. The “spacewatch” telescope on Kitt Peak finds several close approaches per 
month. [Ref. 4], 
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Many ECAs have come relatively close to the Earth. Hermes, in 1937, passed 
within 800,000 km (twice the distance to the Moon) of the Earth. It’s orbit was lost and 
it’s location is no longer known. On January 18, 1991, 1991BA passed within 
approximately 170,000 km (0.0011 AU) of the Earth. The newly discovered near-Earth 
asteroid 1994 XM1 made the closest approach to the Earth of any object discovered 
outside the earth's atmosphere—some 105,000 km on the morning of Dec. 9 over Russia. 
The diameter was estimated to be approximately 6-13 meters. 

There are many more asteroids with the potential for inner planet encounters. 
Among them are 433 Eros, 887 Alinda, 1036 Ganymed, 1221 Amor, 1566 Icarus, 1580 
Betulia, 1620 Geographos, and 1627 Ivar. There are many other inner asteroids of note. 
3200 Phaethon has the greatest asteroid eccentricity, and the smallest perihelion, 0.135 
AU, 1951 Lick is located entirely between the orbits of Earth and Mars; it doesn’t cross 
any orbit. 2062 Aten is entirely inside Earth’s orbit. 2102 Tantalus has the highest 
inclination of any numbered asteroid (64°). 1973na has an even greater inclination (68°). 

4. Earth-Crossing Asteroids 

There are many asteroids well within the orbits of the inner planets. A current list 
of all ECAs through 1991 in contained in reference 5. These asteroids are of interest 
because they lie well inside the main asteroid belt. Some of these come quite close to the 
orbit, of Earth. An Earth-crossing Asteroid (ECA) has been defined as a minor planet 
whose orbit can intersect the capture cross-section of the Earth as a result of secular 
gravitational perturbations. [Ref. 5] ECAs have secular periods on the order of tens of 
thousands of years. 

a. Secular Variations 

It has been shown that an asteroid orbit which overlaps the orbit of a planet 
may, as a result of the advance of the line of apsides, intersect the orbit of the planet. If 



the orbit of the planet is of low eccentricity and is overlapped by an asteroid orbit which is 
highly eccentric, the two orbits will be linked, like links of a chain, when the argument of 
perihelion is 0 or 7t and unlinked when the unlinked at %I2 or 3n/2. In one complete 
rotation of the line of apsides, the orbits transition from linked to unlinked states a total of 
four times, providing four possible intersections of the two orbits. These intersections of 
crossings may be found by solving the polar equation of the ellipse, and occur at 



0 ) 


where o is the argument of perihelion, a is the semi-major axis of the asteroid orbit, e is 
the eccentricity of asteroid orbit at the time of intersection, and p is the radius to the 
planet’s orbit along the line of nodes at the time of intersection. [Ref. 6] Because most 
planet’s orbits are not circular and there are secular variations in e, p will have four 
different values for each rotation of the line of apsides which can be found by 
simultaneously solving Eqn (1) and the polar equation for the elliptical orbit of the planet. 
This type of orbit is referred to as a quadruple crosser, and is shown in Figure 2 using 
asteroid 2062 Aten as an example. The figure shows one cycle of advance of q starting at 
(0 = 0. The solid line shows the ascending node and the dashed line shows the descending 
node. 



Figure 2. Secular Variation of Radius to the Node of the Orbit of 2062 
Aten. From Ref. 6. 
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Secular variations in e may be sufficiently high to change the depth of 
orbital penetration of the assailant body. In this case the orbits may transition from a non- 
overlapped condition at a = 0, Jt to a relatively deep overlap at an/2, 3n/2. The first 
crossing occurs due to the secular increase in e as a begins its rotation from 0 to rt/2, thus 
linking the orbits. The second crossing occurs as a increases sufficiently to unlink the 
orbits. Thus, there are two crossings during each 72 advance of the line of apsides. This 
type of asteroid is known as an octuple crosser. An example is shown in using asteroid 
1580 Betulia, the first asteroid where this type of behavior was recognized. [Ref. 6] 
Approximately 1.5 cycles of advance of a are shown, with the solid line representing the 
ascending node and the dashed line representing the descending node. 



Figure 3. Secular Variation of Radius to the Node of the Orbit of 
1974MA. From Ref. 6. 

If the inclination of the asteroid with respect to the ecliptic plane is 
sufficiently high, the secular perturbations are not strong enough to influence a 
continuous advance in a. In this case, a librates around nil or 37i/2. When this occurs, a 
large oscillation of e can occur during the libration cycle leading to orbit intersection four 
times during each libration cycle of a. These types of asteroid are called quadruple 
crossing w librators. An example is shown in Figure 4 using asteroid 1973 NA. 
Approximately five libration cycles of omega are shown. [Ref. 6] , 
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Figure 4. Secular Variation of Radius to the Node of the 
Orbit of 1973 NA. From Ref. 6. 


The fourth type of ECA are known as supercrossers and result from 
asteroids that librate about the 3:1 commensurability with Jupiter. The resonant effects of 
these perturbations cause a relatively high frequency oscillation in a and e which lead to 
relatively high frequency oscillations (4 cycles in approximately 1400 years) between 
periods of overlap and nonoverlap of the two orbits. For example, asteroid 1915 
Quetzalcoatl, shown in Figure 5, had nine crossings of the Earth’s orbit over a period of 
approximately 1400 years. Approximately four cycles of libration of the mean motion of 
the asteroid are shown. [Ref. 6], 



Figure 5. Secular Variation of Radius to the Ascending 
Node of 1915 Quetzalcoatl. From Ref. 6. 

b. Classes of Earth-Crossing Asteroids 

In addition to classifying asteroids by the type of secular perturbations they 
experience. Earth-crossing asteroids are further divided into three groups on the basis of 
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their present osculating orbital elements: (1) Aten asteroids, (2) Apollo asteroids, and 
(3) Earth-crossing Amor asteroids. 

(1) Aten. The orbits of Aten asteroids have a semi-major axis less 
than 1 AU. Their orbits overlap those of the Earth near their aphelia, as shown in Figure 
6 . 

View Looking Down +Z 



Figure 6. An Aten Type Asteroid. 

This class includes all asteroids with a < 1.0 AU and with the aphelion distance of the 
asteroid, Q > 0.983 AU. The asteroid’s orbit can intersect that of Earth as the Earth’s 
perihelion distance is 0.983 AU. The first Aten was discovered relatively recently, in 
1976. Fifteen were known as of August 1993. They represent approximately 10% of the 
180 known ECAs. The first three discovered exhibit continuous, or nearly continuous 
orbital overlap with Earth and are characterized as deep quadruple crossers. [Ref. 6], 

The total number of Aten asteroids to visual magnitude has been 
estimated to be on the order of 100. In the same way that some Earth-crossing Amors 
have current osculating orbits entirely outside the orbit of the Earth, some Atens may 
have orbits entirely inside the orbit of the Earth. Shoemaker et. al. numbers these at a few 
tens of objects out to visual magnitude 18. [Ref. 6] 
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(2) Apollo. The orbits of Apollo asteroids are those asteroids with 
a-> 1.0 AU, and perihelion distance, <7 < 1.017 AU, where 1.017 is the aphelion distance of 
the Earth. An example is presented in Figure 7. 

View Looking Down +Z axis 



125 Apollos have been discovered, which represents approximately two-thirds of all 
ECAs. All but about 5% of all Apollo asteroids are ECAs. [Refs. 6 and 5] Of those 
known in 1979, approximately 60% were quadruple crossers, 20% were octuple crossers, 
5% were quadruple crossers part of the time at octuple crossers part of the time, and 15% 
were quadruple crossing co librators. [Ref. 6] In 1932 the first Earth-crossing asteroid, 

1862 Apollo, was found by K. Reinmuth at Heidelberg. [Ref. 8] 

(3) Amor. Although Amor asteroids have perhaps the simplest 
definition of all, the fact that some are in Earth-crossing orbits can be confusing. Their 
orbits are defined strictly upon the orbital perihelion distance, 1.017 < <7 2 < 1.3 AU. They 
have perihelions close, but a little outside of the Earth’s orbit. The first discovery of an 
asteroid of this class was 433 Eros in 1898. [Ref. 4] However, the traditional name for 
Atens comes from the asteroid of the same name discovered by E. Delporte at Uccle, in 
1932. [Ref. 8] The 47 known Amor asteroids comprise approximately 25% of all ECAs. 
[Ref. 5] The upper limit for perihelion distance of 1.3 AU was chosen because is it is 
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near a minimum in the radial frequency distribution of q for discovered objects. [Ref. 6] 
As discussed previously, some Amors can become Apollos and vice versa, due to secular 
perturbations. The classifications are based on the category at the time of discovery. For 
example, if the Amor asteroid 1915 Quetzalcoatl had been discovered in 1942, it would 
have been classified as an Apollo asteroid. [Ref. 6] An example of an Earth-crossing 
Amor is shown Figure 8. 





c. Distant Asteroids 

While most asteroids are contained within the main belt, there are some 
that lie outside that area. The Trojans, shown in Figure 1, are a group of approximately 
100 asteroids in two thinly populated lobes centered on Jupiter’s orbit at about the 
Lagrangian L4 and L5 points. They are in stable orbits and have an orbital resonance of 
1:1 (the ratio of the number of orbits of the asteroid to that of a third body) with Jupiter. 
While not tightly grouped about the Lagrangian points, they librate about these nominal 
positions, moving closer and further from Jupiter, but not far from its orbit. [Ref. 4] 
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Opposite Jupiter, there is a thinner lobe centered on the Lagrangian L3 
point. These are known as the Hilda asteroids and have a 3:2 resonance with Jupiter. 

Over 40 Hilda type asteroids have been discovered. Hilda orbits are very similar to some 
short-period comet orbits. The difference being that the Hilda asteroids are in stable 
orbits while the short-period comets are not. [Ref. 4] 

Another asteroid of note is 279 Thule. It is the only known 4:3 asteroid. It 
is unique in that having a resonance near unity places it near Jupiter’s orbit where it is 
strongly perturbed every three Jovian years, yet remains in a stable orbit. 944 Hidalgo has 
a comet-like orbit that takes it out nearly to Saturn. 1373 Cincinnati experiences 
relatively strong perturbations from Jupiter but has no resonance making it unique. 

2060 Chiron was the first trans-Saturn asteroid. Discovered in 1977, it has 
a diameter of approximately 200 km. Its motion has been simulated by researchers back 
to 1664 BC, when it appears to have come within approximately 0.1 AU of Saturn. 
1992QB1 appears to be of similar size to Chiron but its orbit extends beyond that of 
Pluto. Exotic bodies such as these may be from the long-postulated Kuiper belt, 
discussed in the comet section. 

d. Asteroid Classifications 

Asteroids are of four different classifications based upon the photometric 
characteristics of the reflected light, the amount of which is a function of their nature and 
composition of their surface material. This classification is only based upon spectral 
characteristics. The main types have good correlation to location in the asteroid belt. 

Most asteroids fall into one of the following four categories. [Ref. 4], 

S-type. A broad distribution of asteroids centered at about 2.3 AU 
from the Sun. This type is associated with stony-iron meteorites, although some scientists 
question this association. These asteroids have high concentrations of nickel and iron, 
with approximately equal amounts of metals and silicates. These are presumed to be the 
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exposed fragments from the inner cores of larger parent bodies. These are the second 
most prevalent type of asteroid. [Ref. 4]. 

C-type. The most common type of asteroid with a relatively sharp 
distribution centered at about 3.1 AU from the Sun. Thought to be composed of material 
similar to carbonaceous chondritic meteorites, C-types are very dark with an albedo of 
approximately 0.04. They have various grainy, carbon rich, rock-like mineralogies, and 
are a complicated group to examine photospectrally. Most scientists agree Ceres is a C- 
type asteroid. [Ref. 4]. 

M-type. Asteroids with moderate reflectivity and a reddish spectra 
associated with metals, particularly Ni-Fe are known as M-type asteroids. Similar 
meteorites are the nickel-irons and the enstatite chondrites, which have Ni-Fe embedded 
in silicates. There are not many M-type asteroids. The distinction between M-type and 
S-type asteroids is not clear and there are many similarities between them. [Ref. 4], 

D-type. D-type asteroids have low albedos and reddish spectra and 
are primarily the Trojans, located near the orbit of Jupiter. They are composed mostly of 
clays with magnetite and carbon-rich materials. No meteorites match these spectral 
characteristics but the dark side of Saturn’s Iapetus is a good match. [Ref. 4]. 

e. Earth-Crossing Asteroid Population 

Visual magnitude is defined as an arbitrary number, measured on a 
logarithmic scale, used to indicate the brightness of an object. A one magnitude 
difference is the fifth root of 100, and is approximately equal to a factor of 2.512. The 
brighter the asteroid, the lower the numerical value of magnitude. Very bright objects 
have negative magnitudes. 

The absolute magnitude, H, of an asteroid is the magnitude the asteroid 
would have if it were 1 AU from the Sun, and viewed from the Sun. For example, if an 
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asteroid was 2 AU from the Sun and viewed from 1 AU away, it would be four times 
dimmer, or about 1.5 magnitude less bright. 

As most ECAs are discovered through direct telescopic observation, the 
absolute magnitude of the asteroids is of extreme importance to an observer. In general, 
the larger the diameter of an object, the lower the value of absolute magnitude (the object 
is brighter than a similar asteroid of smaller diameter). As a result, the large asteroids are 
often discovered more quickly than the smaller objects. 

The two most common types of asteroids are C-type and S-type. Table 2 
presents the current estimates of EC A discovery completeness based on absolute 
magnitude, H. Estimated diameters, based on the albedo of the two most prevalent types 
of asteroids, are also presented. 


H 

Percent 

discovered 

Asteroid type 

C-type diameter 

S-type diameter 

13.2 

100% 

12 km 

6 km 

15.0 

35% 

6 km 

3 km 

16.0 

15% 

4 km 

2km 

17.7 

7% 

2 km 

1 km 


Table 2. Estimated Discovery Completeness From Ref. 5 


Numerous studies have attempted to estimate the total number of asteroids 
in Earth-crossing orbits. These estimates are based on a power law relationship 

N = k D b (2) 


where N is the number of asteroids larger than a given diameter, D. £ is a constant and b 
is the power-law exponent. The general form of the size distribution is based on 
observation as the actual detailed distributions remain unknown. 

There are approximately 180 known Earth-crossing asteroids. The 
characteristics of Apollo asteroids that overlap the Earth’s orbit part of the time and the 
characteristics of Earth-crossing Amors are essentially identical. The distinction is 


19 






derived from the present state of the cycle of variation of their perihelion distance. The 
majority of Earth-crossing Amors are shallow crossers, while the majority of Apollos are 
deep crossers. Because the orbits of most Earth-crossing Amors only overlap the orbit of 
Earth for a small period of time, they retain their traditional classification as Amor 
asteroids. [Ref. 6] 

Based on the best information to date, an estimate of the number of Earth¬ 
crossing asteroids larger than a given diameter, D, are shown in Figure 9. For this 
population model, changes in the power law are estimated to occur at diameters of 10 m, 
70 m, and 3.5 km. 



Figure 9. Estimated Number of Earth-Crossing Asteroids. From Ref. 5. 

C. COMETS 

Comets are a diffuse bodies of gas and solid particles (such as CN, C2, NH3, and 
OH), which orbit the Sun. Their orbits are highly elliptical or even parabolic in nature. 
Edmund Halley used Newton’s celestial mechanics to show that comets orbited the Sun. 
He then deduced that one spectacular comet in particular was returning to Earth with a 
period of 76 years. His correct prediction of the comet’s return in 1758 made it comet 
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Halley. Today, as many as 10 n to 10 12 comets remain, orbiting the solar system in the 
Oort Cloud and Kuiper Belt, described in Section 2, below. 

1. Comet Classifications 

Comets which have an elliptical orbit are known as periodic comets. A significant 
number of comets also have parabolic or hyperbolic orbits and will, therefore, make only 
one perihelion passage. The periodic comets are further divided into short-period and 
long-period comets. 

a. Short-Period Comets 

Short period comets are defined as those whose orbits lie predominately 
within the realm of the solar system. These include all the comets which have been seen 
more than once. The shortest known period is that of comet Encke at 3.3 years while the 
longest extend up to 200 years. Most short-period comets have inclinations close to the 
ecliptic, less than approximately 35°. [Refs. 4 and 3] About 95% move in a direct sense. 
Comet Halley is a short period comet in a retrograde orbit. More recently, short-period 
comets have been further divided into two additional classes: (1) Jupiter and (2) Halley. 
Comets in the Jupiter family have periods of less than 20 years with inclinations close to 
that of the ecliptic and direct orbits. Halley type comets have period of 20 < P < 200 
years. In general, Halley type comets have higher average inclinations than Jupiter type 
comets. 


b. Long-Period Comets 

Long period comets are those with periods greater than 200 years. Long- 
period comets are randomly distributed on the celestial sphere. Often it is difficult to 
distinguish an orbit with a very long period (on the order of 10 s - 10 7 years) from a 
parabolic or hyperbolic orbit. Many appear to have aphelia of 20,000-50,000 AU. These 
types of comets come from all directions. There is no relationship for the orientation of 
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the orbit with respect to the ecliptic or for direct revolution. Approximately one-third of 
all long-period comets observed are on weakly hyperbolic orbits. 

2. The Origin and Nature of Comets 

Many scientists believe comets to be planetesimals formed in the colder reaches 
of the solar system, beyond Saturn, and perhaps beyond Neptune. Those “cometesimals” 
in the planet realm of the solar system were presumed to be ejected by perturbations from 
the outer planets. These comets now reside in the Oort Cloud, at the outer reaches of our 
Sun’s gravitational influence. The Oort cloud, named after Jan Oort in 1950, is a 
reservoir of comets located about a light-year from the Sun. 

Comets are brought into the Earth’s realm from their orbits in the Oort cloud or 
the Kuiper belt as they are perturbed by Jupiter or passing stars every few million years 
and result in comet insertion back into the solar system. [Ref. 9] 

At the same time Jan Oort was formulating his hypothesis, the structure of comets 
was proposed by Fred Whipple. His “dirty snowball” model hypothesizes a composition 
of various ices, predominantly H2O, along with lesser amounts of silicate and other 
mineral dusts. This analysis of the dirty snowball/Oort cloud has remained the generally 
accepted model, with subsequent modification and elaboration due to follow-on research. 
Relatively recent comet space probes, particularly to comet Halley in 1986, have greatly 
increased general knowledge about the comet. 

3. The Oort Cloud 

The factors that led to Oort’s comet cloud model were 1) very long-period orbits 
that were commonly found to have aphelia of around 50,000 AU, 2) the large reservoir 
that would be required in order to supply comets over the last 4 billion years at the rate 
they appear in the inner solar system (a trillion or so), and 3) the fact that very long- 
period comets appear from all directions. Oort proposed that the planetary system was 
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surrounded by a distant spherical cloud of comets extending approximately one-half the 
distance to the nearest stars. The comet cloud is perturbed by passing stars which cause 
the injection of a long-period comet into our solar system. The source of the Oort cloud 
is postulated to be icy planetesimals ejected by the proto-planets in the outer solar system, 
particularly Uranus and Neptune. [Ref. 3] 

4. The Kuiper Belt 

The Kuiper Belt is a disk thought to be composed icy planetesimal remnants 
beyond Neptune, and was first proposed by Kuiper in 1951. Because these remnants have 
such a high orbital periods and because the density of material in the solar nebula 
accretion disk decreases beyond Neptune, this material was unable to accrete into a 
planetary body. It has been postulated that the Kuiper Belt is responsible for supplying 
the low-inclination Jupiter-class comets. [Ref. 3] 
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III. DEFENDING THE EARTH AGAINST IMPACTS FROM 
LARGE COMETS AND ASTEROIDS. 


A. INTRODUCTION 

There are numerous illustrations of asteroid and comet impacts on Earth: 
meteorite craters, astroblemns, and certain circular geologic structures. Many hypotheses 
have been formulated to predict the consequences of a collision with a Near Earth Object 
(NEO) with the Earth. Damage could be localized to a small geographical area, result in 
regional political instabilities, or cause devastating climatic change, dwarfing the feared 
“nuclear winter” damage of a global thermonuclear war. Recent research into the effects 
of a collision have fairly well established the evidence for severe climate alterations. The 
best known and most extensively investigated impact occurred on the Yucatan Peninsula, 
at the boundary of the Carboniferous and Paleogeneous Ages, 65 million years ago, and is 
thought to have caused the extinction of giant reptiles of that time. The geochemical and 
paleontological record has demonstrated that a 10- to 15-km diameter NEO impacted the 
Earth with the force of 100 million megatons in explosive energy. [Ref. 1] The 
frequency of such impacts is small but not insignificant. 

The most energetic event of this century occurred over the Tunguska Valley in 
1908. It was surmised that this event resulted from the break-up of a chondritic asteroid 
at an altitude of several kilometers. The explosive yield was originally estimated at 15— 
20 megatons (Mt) of TNT, although recent estimates have placed that number as high as 
48 Mt, which leveled approximately 2,000 square kilometers of forest. [Ref. 9]. 

There are numerous ways of deflecting an Earth-crossing asteroid or comet. One 
involves a direct, kinetic energy type weapon, which imparts a change in momentum to 
the assailant object. A second way of mitigating the threat from an ECA or comet is 
through the use of nuclear explosive devices. Issues of great importance surrounding the 
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use of a nuclear device on a potential colliding object include: 1) estimating the 
radiation energy transfer to its surface; 2) maximizing radiation energy transfer; 3) 
understanding the hydrodynamic motion and mass blow-off from the surface, and the 
resultant momentum transfer; and 4) gaining knowledge of the various ways the object 
may be disrupted due to fracture or fragmentation. [Ref. 9]. Regardless of whether one 
of these, or some other method, is chosen, the several basic issues must be thoroughly 
understood before an attempt is undertaken to deflect an object. 

The most likely candidate for mitigating the danger of an Earth impact is a rocket- 
delivered nuclear explosive. Nuclear explosive devices possess the highest concentration 
of energy, and can be manufactured with warheads yielding 100 Mt or more. A 
thermonuclear charge with a yield of 1 Mt has a mass of the order of 0.5 ton. Of interest 
is the fact that to obtain the equivalent energy by impact of a body of the same mass 
requires an impact velocity of approximately 4,000 km/s. [Ref. 9]. 

Simonenko et al. [Ref. 9] presents a list of problems that must be considered 
before the development of a defensive system could be undertaken. 

1. Assessment of the number of dangerous cosmic bodies and their probability of 
their collision with the Earth. 

2. Assessment of the time required to detect and identify them as dangerous. 

3. Understanding of the object’s motion in the immediate vicinity of the Earth, 
penetration of the atmosphere, and detailed dynamics during collision with the Earth, 
with concomitant assessment of the local and global consequences of the collision. 

4. Assessment of the required effect on the astral assailant, whether it be 
deflection or fragmentation. 

5. Consideration of the needed nuclear devices, delivery means, and optimal 
regime of action. 

6. Assessment of consequences of collision with fragments of an object that had 
been fractured by a nuclear explosive. 
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7. Consideration of ecological consequences for the Earth and space environment 
of nuclear explosions in space. 

The most important question from a technical viewpoint is to ascertain if it is 
feasible to significantly alter the trajectory of an ECA or comet by detonating a nuclear 
explosive near the surface. Equally as important is the ability to calculate the effects of 
such a weapon on the extraterrestrial body. Many characteristics of nuclear explosions 
are predictable and well defined. For example, detonation of a 100-kt device known as 
“Sedan” produces a crater 370 m in diameter and 100 m in depth. Significant differences 
may result by varying the nuclear explosion yield, distance to the target object, the object 
dimensions, material composition of the object, and the design of the nuclear device. 
Depending upon the interplay of these variables, the resultant explosion can result in 
fragmentation, crushing, or deviation of the target object from its initial trajectory. [Ref. 

91 - 

Table 3 presents the yield vs. mass for nuclear explosive devices. This shows that 
the yield from current nuclear technology, that which can be launched with existing 
rocket technology, is capable of deflecting small objects. The nuclear device can be 
scaled upwards, perhaps an order of magnitude or more, while preserving special 
characteristics. In this case, modification would not require further testing. [Ref. 9] Of 
particular importance to the deflection of an asteroid are the specific physical 
characteristics of the body. These include shape, material composition, and physical 
stability of the object. Maximum effectiveness for deflection may require several nuclear 
devices detonated in a carefully determined sequence, or one device of special 
configuration. In addition, special measures can be taken to minimize the radiation 
contamination of the asteroid fragments, particularly important if it is determined that the 
fragments from such an explosion will continue towards impact with Earth. 
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Yield 


Mass 


1 Mt 
10 Mt 
100 Mt 


0.5 ton 
3 to 4 ton 
20 to 25 ton 


Table 3. Yield vs. Mass for Nuclear Devices [Ref. 9] 


Finally, an area of study which needs additional study is that of the explosion of 
nuclear devices in space. Major differences include absence of atmosphere, 
commensurability of the object’s dimensions with linear scales of the phenomenon, 
complexity of object shape, relatively weak gravity, and exotic material composition. 
[Ref. 9] 

B. DEFLECTION CONSIDERATIONS 

There are two main areas of interest in the mitigation of a Near-Earth Object 
(NEO): 1) causing the reatening object to break up, and 2) altering the trajectory of 
the object. There are several issues surrounding these areas which concern the scientific 
community with regard to deflecting or fragmenting a NEO sufficiently to avert collision 
with Earth. First, is there a limitation to the size of asteroid that can be deflected, given 
today’s technology. Secondly, given the option, what is the most effective way to deflect 
these astral assailants to ensure, with a reasonable degree of certainty, that they will miss 
the Earth. 


1. Size Limitations 

Simonenko et al. [Ref. 9] has estimated the maximum dimensions of an object 
that can be influenced by a nuclear explosion as follows. Assume a nuclear explosion 
occurs on the surface of an asteroid having mass m a and radius R a . Assume the mass m 
of the material ejected by the explosion is small compared with the mass of the object. 
The mass of material ejected from the body has a distribution of velocities, but to 
simplify the estimation, it is assumed that all the ejected material has the same velocity, 
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the average of all the actual velocities, v m . This varies depending upon the explosive 
yield and material properties of the object. 

Assuming m«m a , conservation of momentum shows 

mv,„ = m a Av a (3) 

where v» is the velocity of the ejected material an infinity and Av o is the velocity increase 
of the object as a result of the explosion. In similar fashion, conservation of energy 
shows 



where G is the universal gravitational constant. Combining Eqns. (3) and (4) gives 



Examination of Eqn. (5) shows that if an asteroid is of sufficient mass, and given 
the above limitations on ejected material velocity, ejected material will simply fall back 
to the asteroid or remain in orbit around it, resulting in the object resuming its original 
velocity. To approximate the maximum radius and mass that could be deflected, Eqn. (5) 
solves for R aC rit with v« = 0. At this radius, all attempts to deflect the object are 
impossible by this method of explosion. The critical radius is defined as 


Ken, 



( 6 ) 


where p a is the average density of the object and 



(7) 


Finally, the critical mass is given by 


2G J 8np a G ' 


( 8 ) 
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Knowles and Brode have published data showing that a conservative estimate of the 
ejection velocity, using currently available technology, is approximately 100 m/s. Using 
and average density of 5 g cm' 3 for the assailant object, an average ejection velocity of 
100 m s' 1 for current technology, the critical mass and radius for the object are 
approximately 5 x 10 18 kg and 60 km, respectively. Asteroids with masses or radii 
greater than these values cannot be deflected with current technology, regardless of the 
warning time involved. [Ref. 9] 

2. Directional Considerations 

The orbit perturbation requirements to deflect objects from collision with the 
Earth are an important consideration when planning a deflection mission. Two scenarios 
are considered: (1) A short time-scale deflection where equations of motion are 
linearized, and (2) A long time-scale deflection, using orbital equations of motion. These 
two scenarios are developed fully below, the results of which will be used to generate the 
change in velocity requirements for a hypothetical asteroid. 

a. Short Time-Scale Deflection 

The short time-scale deflection scenario is easily developed using 
rectilinear equations of motion. This is valid if the time scale involved is short compared 
to the orbital period, P. To prevent a collision, a change in velocity must be imparted 
orthogonal to the flight path of the object. The displacement, 8, that must be achieved by 
a change in velocity, Av, is 

8 = Av • /. (9) 

Assuming the object were on a path such as to intersect the center of mass 
of the Earth, the object must be deflected a minimum of one Earth radius. This provides 
no margin of safety, however, and therefore, taking the minimum acceptable distance to 
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be two Earth radii, the minimum change in velocity required to deflect a potential threat 
is given by 

2 IL 147.6ms -1 (10) 

Av «-a---. 

t f,days 


b. Long Time-Scale Deflection 

For a deflection to be accomplished over a longer period of time, the 
equations of motion cannot be simplified to the rectilinear case and the two-body orbital 
mechanics equations of motion apply. There are three directions with respect to the flight 
path which will cause an orbiting body to deviate from the original orbit: (1) orthogonal 
to the flight path, along the line parallel to the angular momentum vector, (2) orthogonal 

to the flight path, in the plane of the orbit, and (3) along the flight path, in the orbit plane. 

A two body circular orbit is assumed to simplify the problem. 

Given a longer response time, the best deflection direction may not be 
orthogonal to the flight path direction, as in the short time-scale deflection. Each option 
is assessed below to determine the optimal deflection direction. 

(1) Out of plane, orthogonal to flight path. A change in velocity 

made orthogonal to the orbit plane, as shown in Figure 10, 

results in a change of inclination according to 

. „ • * (11) 

Av = 2v sin —, 
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where v is the orbital speed of the object and i is the required inclination change. The 
maximum displacements occur n/2 and 3 ji/ 2 around the orbit. The orbital period remains 
unchanged. Using a small angle approximation gives 


where 8 is the amount the body’s orbit will be deflected n/2 past the point of the 
deflection, and r is the radius of the object’s orbit around the Sun. Finally, substituting in 
for tire orbi tal velocity in terms of orbital period, P, gives 



The minimum change in velocity required to deflect the asteroid two Earth radii is given 
by 
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( 15 ) 


For a circular orbit at 1 AU, the minimum change in velocity required to achieve the 
required deflection is approximately 2.54 m/s and should be imparted 90° prior to the 
expected impact point in the asteroid’s orbit. 

(2) In plane, orthogonal to flight path. The change in velocity 
required to deflect an object using a deflection in the asteroid’s orbital plane, orthogonal 
to the flight path is found in a manner similar to (1), above. As before, the orbital period 
remains the same. The asteroid is displaced along it’s track with the maximum 
displacement occurring 7t/2 and 37t/2 around the orbit, as shown in Figure 11. The 
maximum deflection 



Figure 11. Deflection In Plane, Orthogonal to Flight Path. 


g _ 2AvP ( 16 ) 

max ~ n 

and the minimum change in velocity required to deflect the asteroid two Earth radii is 
given by 
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Av s 


The required change in velocity for this case is approximately 0.635 m s" 1 . 

(3) In plane, along flight path. A deflection along the flight path, 

shown in 



will result in a net increase or decrease in the orbital speed. If there is a net increase in 
speed, the point at which the deflection occurred will be the perihelion of a new orbit. 
Conversely, if there is a net decrease in speed, the deflection point will be the aphelion of 
a new orbit. Irrespective of deflection direction, the period of the deflected object’s orbit 
will change. This change in period serves as the mechanism by which, over time, an 
object can be influenced so as to miss a collision with the Earth. The amount the asteroid 
is displaced during each orbit is directly proportional to the change in period of effected 
by the deflection. This is given by 

S^„=vAt. (18) 

The approximate displacement for a deflection along the flight path can be derived from 
the equations for the specific mechanical energy of a body. 



the orbital period of a body. 
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( 20 ) 


r=2 *llf- 


and the circular velocity for an orbit, 



Taking the derivative of Eqn (19) and substituting v 2 for (i/a gives 
Av Aa 
~ = 2a’ 


( 21 ) 


( 22 ) 


while taking the derivative of Eqn. (20) and resubstituting in for T gives 
A a 

AT = 37—. 


Substituting (22) and (23) into (18) gives the maximum deflection per orbit, in terms of 
the original orbital period of the asteroid and the change in velocity imparted to the 
object, as 

S per _ orbll = 37Av . (24) 


So, the change in velocity required to deflect an object two Earth radii is approximately 


(25) 


where n is the number of orbits prior to impact. If n = 1, the change in velocity is 
approximately 0.135 m s' 1 . This is better than either case (1) or (2). If the impact can be 
predicted as much as a decade in advance, the change in velocity required to deflect the 
asteroid is reduced to 1.35 cm s' 1 . This is an order of magnitude less than case (2) and 
two orders of magnitude less than case (1) and is an achievable quantity in terms of 
technological feasibility. 
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IV. DEVELOPMENT OF HYPOTHETICAL ORBITS 


A. INTRODUCTION 

In order to accurately test the sensitivity of the change in velocity requirements 
presented in section III, accurate orbital parameters for hypothetical asteroid objects were 
first developed. This section first describes the coordinate systems used and the 
characteristics of circular, elliptical, parabolic, and hyperbolic orbits. Next, this section 
describes the process by which the Earth’s position in heliocentric coordinates, at any 
given time, are used to generate the heliocentric coordinates of an intersecting asteroid 
orbit. The size, shape, and inclination of the asteroid orbit are variables defined by the 
user with the only limitation being that the two orbits intersect. MATLAB scripts which 
generate these procedures are contained in the Appendix. Orbits can be generated for 
circular, elliptical, parabolic, or hyperbolic trajectories. 

Development of hypothetical orbits followed a rigorous routine allowing for 
repeatability of the work done while minimizing the workload for completing multiple 
simulations. The methodology and theory used in the development of the MATLAB 
scripts are discussed below. 

B. COORDINATE SYSTEMS 

The first requirement for describing an orbit is to define a suitable reference 
frame. In many cases this means finding an appropriate inertial coordinate system. For 
orbits around the Sun such as planets, asteroids, and comets, the heliocentric-ecliptic 
coordinate system is often used. Satellites in orbit around the Earth use the geocentric- 
equatorial system. [Ref. 10] Another convenient coordinate frame used for describing the 
orbit of a body is the perifocal coordinate system. Each rectangular coordinate frame is 
defined by specifying the origin, the fundamental plane (i.e. the X-Y plane), and the 
direction of each axis. In developing the hypothetical orbits used for the simulations, the 
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heliocentric-ecliptic and perifocal coordinate systems were used extensively. These are 
described in detail below. 

1. Heliocentric Ecliptic Coordinate System 

The heliocentric-ecliptic coordinate system, as the name implies, has its origin at 
the center of the Sun. The letters XYZ describe the three principle axes as shown in 
Figure 13. The fundamental plane, the X-Y plane, coincides with the ecliptic, which is 
the 



Figure 13. The Heliocentric-Ecliptic Coordinate System. From Ref. 10. 

plane defined by the Earth’s orbit around the Sun. The line of intersection of the Earth’s 
equatorial plane and the ecliptic plane (XY) where the Sun crosses the equator from south 
to north in its apparent annual motion along the ecliptic is known as the vernal equinox. 
When the vector from the Earth to the Sun (X-axis) points towards what is now 
approximately the constellation of Pisces, it coincides with the vemal equinox, also 
known as the first day of spring. The Y axis is 90° from the Z axis in the direction of the 
Earth’s rotation around the Sun. The Z axis completes the right-handed coordinate 
system. [Ref. 10] 

2. Perifocal Coordinate System 

The perifocal coordinate system, shown in. Figure 14, is one of the most useful 
coordinate systems for describing the motion of a body. The letters PQW are used to 
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describe the three principle axes. The fundamental plane, the PQ plane, is in the plane of 
the body’s orbit and the origin is at the focus of the primary body, in our case the Sun. 
The P axis is in the direction of perihelion while the Q axis is 90° from the P axis in the 
direction of rotation. The W axis completes the right-handed coordinate system. 



Figure 14. Perifocal Coordinate System. 


C. CLASSICAL ORBITAL ELEMENTS 

Five classical orbital elements, also known as Keplerian elements, are sufficient to 
describe the size, shape, and orientation of an orbit. A sixth element is required to locate 
a body in that orbit. 

There are many ways of locating the position of a body in the orbit plane which 
may be substituted for time of perihelion passage. 

1. True anomaly at epoch, v 0 , is the angle in the plane of the body’s motion 
measured from perihelion to the position to the body at a given time (epoch). 

2. Argument of latitude at epoch, Uo, is the angle in the plane of the orbit between 
the longitude of the ascending node and the radius vector to the body at a given time. If <a 
and v 0 are both defined then 

Uo=co+v 0 . (26) 
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If there is no ascending node, as in an equatorial orbit, then both m and Uo are 


undefined. 

3. True longitude at epoch, 1 0 , is the angle measured between X and r, the radius 
vector to the body, measured eastward to the ascending node, if it exists, and then in the 
orbital plane to r. If Q, co, and v 0 are all defined, then 

lo=f2+co+v 0 =n+v 0 =Q+Uo. (27) 

Two other terms used in the descriptions of orbits are direct, or prograde, and 
retrograde. Direct means easterly, in accordance with the right-hand rule, and is the same 
direction in which the Sun, Earth, and most of the planets and their moons rotate on their 
axes and the direction in which all planets rotate around the Sun. A retrograde orbit, as 
the term implies, is the opposite of direct. Table 4 provides a summary of the various 
types of orbits and the orbital parameters that result. 


Parameter 

Circle 

Ellipse 

Parabola Hyperbola 

Eccentricity 

e=0 

0<e<l 

e=l 

e>l 

Energy 

-|j/2a<0 

-|i/2a<0 

E=0 

mu/2a>0 

Period 

P=27i/n 

P=27t/n 



Anomaly 

Undefined or 
arbitrary 

Eccentric: 

Parabolic, D 
v D 

tan— = , 

2 V2 q 

Hyperbolic, F 


Table 4. Keplerian Orbits 


D. DERIVATION OF FORMULAE 
1. Procedure 

a. Orbital Elements From r and v 

Conversion from position and velocity vectors to orbital elements is done 
using standard procedures. For the purposes of this example, the Earth’s orbital 
parameters are being determined from position and velocity vectors. 
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st the angular momentum is determined. From the relative form of the 


basic two-body equation 


Cross multiplying both sides by the position vector, r, gives 


Expanding this equation, one then finds the angular momentum 


The vector pointing to the node is then calculated using 


where K is the unit vector in the Z direction. If the magnitude of the node vector is zero, 
the orbit is in the fundamental plane. The MATLAB script used to calculate this vector 
has a filter to set the vector equal to zero if the magnitude is less than lxl O' 5 . 

Next the eccentricity vector is calculated with 

The magnitude of this vector provides the eccentricity. It is important to realize that this 
vector is zero for circular orbits. « 

The semi-latus rectum, p, is calculated 


and the radius at aphelion and perihelion 



( 35 ) 


Once the size and shape parameters have been determined, all that remains 
is to find the orientation of the Earth’s orbital plane and the position of the Earth in that 
orbit. First, the inclination is found with 


.’.here is no need to check for angle ambiguities as the inverse cosine always returns an 
angle between 0° and 180°. 

Next, the longitude of the ascending node, which varies from 0° to 360° 
must be determined. By definition, the Earth is in the mean ecliptic and the inclination is 
zero deg. By convention, the longitude of ascending node is defined to be zero deg; there 
is no real “ascending node”. However, due to small perturbations in the Earth’s orbit 
around the Sim, the center of mass of the Earth deviates out of the mean ecliptic to what 
is known as the true ecliptic. As a result, when there are instantaneous components of the 
position and velocity vector normal to the mean ecliptic, very small, but calculable values 
for inclination can be found. Also, a value for the longitude of ascending node can be 
found. If the instantaneous position and velocity vectors lie in the plane of the ecliptic, 
the longitude of ascending node and inclination are zero deg. and the only rotation is that 
of the argument of perihelion. As mentioned previously, orbits aligned with the 
fundamental plane, in our case the ecliptic, have no node vector and the longitude of the 
ascending node is undefined. However, a value for the majority of orbits is obtained from 


cos(Q) = 


M"f 


(37) 


In this case a quadrant check must be conducted and if the dot product of J and n is less 
than zero, then Q must lie between 180° and 360°. 
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The argument of perihelion, which also varies from 0° to 360° is 
calculated. Again, this value is undefined for circular orbits and orbits in the ecliptic 
because perihelion does not exist for the former and there is no node vector in the later. 


~ \n\\e | • 


( 38 ) 


Angle ambiguities are resolved by computing the dot product of K and e. If this value is 
less than zero, then perihelion is below the ecliptic and the argument of perihelion is 
between 180° and 360°. 

The final classical orbital element is the true anomaly at epoch, which can 
vary from 0° to 360°. As mentioned previously, the true anomaly is the angle in the orbit 
plane measured from perihelion to the Earth’s current position. The true anomaly is then 
, , e-r (39) 

C “ (V) ‘W 


The correct quadrant for the true anomaly is found from the dot products of r and v. If 
the dot product is less than zero, then the true anomaly is between 180° and 360°. 


b. Singularity Solutions 

Several special cases must be considered, especially when calculating 
values for the orbital elements using a computer routine. The first occurs when the 
orbital plane is the same as the fundamental plane. This is especially relevant to the orbit 
of the Earth around the Sun. By definition, the Earth is in the ecliptic, the inclination and 
longitude of ascending node are zero and the normal equation to calculate the argument 
of perihelion cannot be used. However, the instantaneous position and velocity vectors 
may contain small components out of the plane of the mean ecliptic. There are small 
values for inclination and longitude of the ascending node. To obtain the required 
accuracy, these angles are included when calculating the asteroid’s orbit. Appropriate 
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filters handle the case when the Earth’s position lies exactly in the mean ecliptic, and the 
inclination and longitude of ascending node are zero. 

These singularities are removed by using the true longitude of perihelion, 
ramie, which combines Q and co. It is found through 




(40) 


wnjvh is valid for elliptical orbits which lie in the ecliptic. Again, a quadrant check is 
f t* i.rsary and is found using the dot product of J and e. If this value is less than zero, 
the true longituue of perihelion is between 180° and 360°. 


c. Orbital Anomaly Determination 

There are several other useful formulas that must be derived prior to 
writing script files for the plotting of orbits and calculation of orbital parameters. First, 
hi elliptical orbits, such as those of the Earth, ECAs, and short-period comets, formulas 
tor determining the eccentric anomaly and the true anomaly prove quite useful. Figure 15 
shows an x,y Cartesian coordinate system centered at C. Fis the focus of an ellipse, in 
this case the Sun. An auxiliary circle of radius a with center C is constructed. P is the 
location of the Earth in its orbit about the Sun and Q is the point where the perpendicular 
to the semi-major axis intersects the auxiliary circle. The angle E is the eccentric 
anomaly and v is the true anomaly of the point P. 
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Figure 15. Orbital Anomalies for Elliptic Motion 


The x and y coordinates of the ellipse can be expressed in parametric form as 

x = acosE y = bsinE. (41) 

for the coordinate system centered at C, shown above. If the focus is selected as the 
center of the coordinate system, the radial position of the point P can be expressed in 
terms of E using the general form of for the radial position of an ellipse with the origin at 
the center 

r = a- ear (42) 

so that 

r = a(l - e cos E). (43) 

Comparing this to the polar equation of the ellipse 

o(l -e 2 ) (44) 

1 + ecosv 


the identities for true anomaly and eccentric anomaly are obtained as 
cos E-e ^ ^ e + cosv 

1 — ecosE 1 + ecosv" 


The true anomaly for the body is obtained by solving Eqn. (44). 

»—-(KM) 


(45) 


( 46 ) 


Another identity. Gauss’ law, provides a useful relation between v and E. 
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tan — 


( 47 ) 



Because 1/2 v and 1/2 E are always in the same quadrant, there is no angle ambiguity to 
resolve when solving this form of the equation. 

2. Determination of r and v From Intersecting Orbit Parameters 

From a given orbit and the position of a body within that orbit, a second orbit can 
be constructed which intersects the first orbit. The location of a body within the second 
orbit which intersects the body in the first orbit is specified by the formulation of the 
solution. Values for semi-major axis, eccentricity, and inclination are variables and can 
be any value with the only limitation that the two orbits must have at least one point of 
intersection. Finally, the position and velocity vectors in the new orbit at the point of 
intersection are calculated. These position and velocity vectors are in perifocal 
coordinates. The detailed procedure for calculating these vectors is presented below. 

a. Circular, Elliptical, and Hyperbolic Orbits. 

The polar equation for an ellipse 

a(l-e 2 ) (48) 

1 + ecosv 

can be used to solve for the true anomaly at epoch. The magnitude of r in this equation is 
the magnitude of the Earth’s position vector at epoch. Once the true anomaly at epoch is 
known, the position vector in the perifocal coordinate system can be determined using 

r = r cos vP + r sin vQ (49) 

where the magnitude of r is once again the magnitude of the Earth’s position vector and 
epoch. 

Differentiating Eqn. (49) in this “inertial” perifocal frame gives 
v = r = (r cosv - rv sinv)P + (r sinv + rv cosv)Q. (50) 
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Substituting ii 


gives 




[u 

l—esinv rv = —(1 + ecosv) 

V P 




sinvP + (e + cosv)Q] 


(51) 


(52) 


b. Parabolic Trajectories. 

The semi-major axis of a parabolic orbit is undefined. If the desired orbit 
is parabolic in shape, the parameter, p, is substituted for the semi-major axis. The polar 
equation of the parabola is then 

r = _P__ (53) 

1 + cosv ' 

This equation is solved for v, as discussed previously. The formulas for 
calculating the position and velocity vectors in perifocal coordinates are then identical to 
the elliptical case discussed previously. 

3. Orbit Rotations 

Rotations between coordinate frames are accomplished using Direction Cosine 
Matrices (DCMs). Coordinate rotations may be accomplished by any number of 
methods. One common method of conducting coordinate rotations is through a 3-1-3 
rotation: A rotation about the 3-axis, followed by a rotation about the new 1-axis, and 
finally a 3-rotation about the new 3-axis. And example, shown in Figure 16, converts 
from perifocal coordinates to Heliocentric-ecliptic coordinates. The rotation matrices as a 
function of angle of rotation, a, are presented below. The subscript associated with each 
transformation is the axis about which the rotation is to be accomplished. 
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C,(a) = 


1 0 0 
0 cosa sina 


( 54 ) 


_0 -sina cosa. 


cosa 


C,(a) = 


-sina 


0 


sina 0 
cosa 0 
0 1 . 


(55) 


Transforming an orbit expressed in perifocal coordinates to heliocentric 
coordinates is conducted through the following formulations: 


= C 3 (co )C, (z')Q (Ci)r re w ■ 

(56) 

= (co )C,(i)C 3 (Q)v pgiy . 

(57) 



Figure 16. Coordinate Transformations. From Ref. 10. 
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V. CASE STUDIES. 


A. INTRODUCTION 

Two tests were conducted to analyze the effects of orbital perturbations on an 
assailant asteroid. The magnitude of the velocity applied, is the change required along the 
flight path, to change the orbital period of a circular orbit such that the change in position 
after one orbit is equal to two Earth radii. This magnitude is approximately 13.5 
cm/sec/orbit. The first used a solar system simulation modeling program called Dance of 
the Planets. [Ref. 4] This program, described in further detail below, allows the user to 
input hypothetical objects into a relatively high fidelity simulation of the solar system. It 
includes gravitational models for the Sun and all of the planets. Objects can be imported 
into the program using their orbital elements or their position and velocity vectors at a 
given time. The model then calculates new orbital positions and velocities for each body 
included in the model using a full gravitational model. The simulation accurately 
computes both forward and backward in time. Although the program provides a 
reasonable model visualization, an inaccuracy was discovered during testing which 
prevented importing of position and velocity vectors of sufficient accuracy for valid tests 
to be conducted. This inaccuracy is demonstrated and explained in section 7 of this 
chapter. The discussion focuses on the procedure used to generate the desired asteroid 
vectors. This would prove useful for follow-on research should the program be modified 
to provide the required accuracy. 

The second method used to evaluate sensitivity of deflection distance and 
direction to variations in orbit eccentricity and size (semi-major axis) consisted of a 
relatively simple model to determine miss distance based on the asteroid’s mean anomaly, 
M, and mean motion, n. The position and velocity vectors at some time prior to collision 
were calculated, new orbital elements determined, and the new position at the original 
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time found. The distance between the center of the Earth’s orbit and the asteroid’s orbit 
was then calculated. The effect of varying specific orbital parameters while keeping the 
magnitude of the change in velocity constant were analyzed. Parameters varied were 
orbital eccentricity, semi-major axis, and deflection direction. 

B. DANCE OF THE PLANETS 

1. Hypothetical Orbit Determination 

The purpose of this section is to present the process used to calculate asteroid 
object orbital elements and heliocentric coordinates from the Earth’s heliocentric 
coordinates. The position and velocity vectors obtained are the position of the asteroid in 
its new orbit around the sun in heliocentric coordinates. This is exactly the same position 
as that of the Earth, hence collision. The asteroids velocity vector is that vector which 
meets the size and shape requirements dictated earlier in the orbit development process 
and is the asteroid’s instantaneous velocity vector at collision. These vectors are then 
imported into Dance of the Planets at the epoch for which the Earth’s original r and v 
were obtained. The orbit was then propagated backwards in time. New velocity vectors 
with a delta-v applied were calculated and imported into the program. These new vectors 
were propagated forward to determine if collision had been averted. The procedure, 
summarized here, is described in detail in the sections which follow. 

a. First, the Earth’s orbital elements are calculated from r and v. 
Specifically, true anomaly, v, is required so that the position in the 
orbit is known. 

b. The asteroid’s orbit size, shape, and inclination (a, e, and i) are 
selected. 

.. The asteroid’s orbital elements are calculated in perifocal coordinates. 

d. The asteroid’s orbit in perifocal coordinates is rotated, using a 3-1-3 
rotation, to make it coplanar with the Earth’s orbit. The eccentricity 
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vectors are aligned. This makes the Sun and the two orbit’s 
perihelions collinear. 

e. A second 3-1-3 rotation is performed to obtain the desired inclination 
and align the two r vectors. The first 3-rotation is through an angle 
Veaith and aligns the asteroid’s perihelion with the position of the Earth 
in its orbit. The 1-rotation inclines the orbit to the desired value. The 
final 3-rotation, through an angle -Vaster, aligns the asteroid with the 
Earth’s current position. 


2. Initialization. 

The first step in the development of hypothetical orbits and completing the 
simulation of Near Earth Objects (NEOs) was to select a desired date and time for the 
impact. For the purposes of this research, the date chosen was January 01, 2007 at 00:00 
UT. This date was selected as it gave a real world flavor to the problem, approximately 
one decade would elapse from the time of discovery to the time of impact. Universal 
Time (UT) is the standard time in the time zone centered on zero degrees longitude in 
Greenwich, England. UT is used for calculating where an object is in the sky relative to 
the horizon at a particular location. 

Next the position and velocity vectors for Earth at the desired impact date were 
obtained from Dance of the Planets. The operator’s manual contains an excellent tutorial 
and provides the information required to operate the program. Several preliminary steps 
must be completed prior to date initialization. The simulation is initialized from Space 
View. Second, from the main access screen, any planets that are to be excluded from the 
gravitational model are deselected. This is indicated by a zero in the on column next to 
the planet name. Returning to the simulation Space View, confirm the simulation pace is 


This minimizes the elapsed time upon returning to the simuk 





initializing the date. In addition, this also minimizes the change in planetary position and 
velocity vectors during the time from initialization to reselecting the main access facilities 
display where the Earth’s vectors are obtained. Finally, enter the desired date. Digital 
hours and minutes appear at low simulation paces. Local time, determined by site 
longitude, is shown unless time is offset. 

Immediately after entering the desired date, reenter the main access screen. For 
each planet selected, Dance of the Planets calculates the period, the distance to the Earth 
and Sun in astronomical units (AU), heliocentric longitude (calculated and simulated), 
right ascension of the ascending node, declination, and apparent magnitude. The 
difference between calculated and simulated heliocentric longitude is that the calculated 
value is the position as determined from ephemeris equations while the simulated value 
gives the position in “simulator space.” Comparison of the two values provides an 
accuracy check for the user between the ephemeris equations and the Dance of the 
Planets model. 

The instantaneous position and velocity vectors of the planets provide the 
essential information to describe dynamic state. With these datum the six orbital 
parameters of eccentricity, semi-major axis, right-ascension of the ascending node, 
inclination, argument of perihelion, and true anomaly, can be calculated. 

3. Constants 

The program provides the user with the option of obtaining the position and 
velocity vectors or the orbital elements at a given epoch. The display shows the 
calculated vectors for the planets in heliocentric coordinates, XYZ, as well as the epoch 
Julian date. Units for the position and velocity vectors are E6 km and km/sec, 
respectively. The epoch date is presented in the form 24ddddd.ffff and is reference to the 
epoch (equinox) of the year 2000 (J2000), as are contemporary star charts. Only Julian 
dates beginning with 24ddddd can be used, giving a 14,680 year window for simulations. 


52 




[Ref. 4] Mean planetary constants for epoch J2000 of the Earth are given in Table 5. 
Solar constants are shown in Table 6. 


Parameter 

Value 

Semi-major axis 

AU 

1.000 001 017 8 

km 

149,598,023 

Eccentricity 

0.016 708 617 

Inclination ~ deg. 

0.000 000 000 

Long, of ascending node (Q) ~ deg 

0.000 000 00 

Long, of perihelion (co) ~ deg 

102.937 348 08 

True longitude (L) ~ deg 

100.466 448 51 

Orbital period (P) ~ years 

0.999 978 62 

Orbital velocity (v) ~ km/s 

29.77859 

Equatorial radius (Re) ~ km 

6378.137 

Gravitational parameter (p) ~ km 3 /s 2 

3.986xl0 5 

Mass (Me) ~ kg 

5.9742xl0 24 

Rotation ~ days 

0.997 269 69 

Inc. of equator to ecliptic ~ deg 

23.45 

Density ~ gm/cm 3 

5.515 


Table 5. Mean Earth Constants for Epoch J2000 


Parameter 

Value 

Radius of the Sun (Rs) ~ km 

696,000.000 

1.0 AU ~ km 

149,597,870.0 

1.0 AU/TUsw, ~ km/solar s 

29.784 691 674 9 

Gravitational parameter (p) ~ km 3 /s 2 

1.327 124 28x10" 

1.0 TU ~ solar days 

58.132 440 9 


Table 6. Solar Constants 


4. Orbital Element Selection 

The orbital elements selected for simulation, shown in Table 7, were selected 
from a list of all known ECA’s given by Rabinowitz et. al. [Ref 5] and [Ref 7], These 
were selected because they each have a point of close approach distance within 20 lunar 


53 






radii over the course of the next century, and because they are representative of each class 
of asteroid as well as a short-period comet. 


No. and Name 


1 Approx, 
diamete 
r(km) 

Depth 

(AU) 

q 

(AU) 

-71- 

(AU) 

— 

-jT- 

(deg) 

2340 Hathor 

20.26 

0.2 

0.356 

0.464 

0.844 

0.450 

5.85 

4660 Nereus 

18.30 

1 

0.092 

0.954 

1.490 

0.360 

1.43 

4179 Toutatis 

14.0 

5 

- 

0.542 

2.505 

0.640 

0.47 

Encke 

9.8 

- 

- 

0.331 

2.283 

0.855 

11.9 


Note: (1) Values used in the development of the hypothetical orbit 


Table 7. Selected Near-Earth Objects Used for Simulation Study. 


5. Asteroid Orbit Rotations 

Two series of 3-1-3 rotations are used to achieve the desired orbit. The first series 
of rotations transforms the position and velocity vectors from perifocal coordinates to the 
Earth’s orbital plane. These angles were previously calculated and are the Q, i, co rotation 
angles. Because the true ecliptic is used, the Earth’s orbit with respect to the mean 
ecliptic is found using the three rotations calculated in Eqns. (37), (38), and (39), 
inclination, longitude of ascending node, and argument of perihelion, respectively. 

Once the assailant object’s orbit is coplanar with that of Earth, and the eccentricity 
vectors are aligned, a second 3-1-3 coordinate transformation is accomplished. This 
series of rotations serves two important purposes: (1) they rotate the orbit to the desired 
inclination, and (2) they align the position vectors of the Earth and the asteroid. 

The first 3-rotation, through the angle v carth , aligns the object’s eccentricity vector, 
which points towards perihelion, with the Earth’s position vector. The 1-rotation, about 
the eccentricity vector, effects the desired value of inclination. The final 3-rotation about 
the orbit normal, though an angle -Vasiami, aligns the asteroid’s position vector with the 
Earth’s position vector. 
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6. Simulatioi 


Once the position and velocity vectors are calculated, they are inserted into Dance 
of the Planets as object vectors. The fictitious orbit, at epoch, is one having the same r as 
Earth but with a velocity vector as provided by MATLAB to provide the desired orbit size 
and shape. 

The simulation is run backwards in time to an arbitrary time. This will he the new 
start time of the simulation. At this point, the simulation is paused, and the new orbital 
elements are calculated for the asteroids current orbital set. These are not the same as the 
original orbit’s parameters due to the influence of third bodies on the asteroid’s orbit 
during the simulation period. Specifically, the influence of the Earth on the asteroid’s 
orbit during the first month of the reverse simulation has a significant impact on the orbital 
parameters. 

Next, another MATLAB script file uses the current orbital parameters to 
recalculate the velocity vector of the assailant body for a desired delta-V applied to the 
body at a given time. This procedure is described in detail below. 

Finally, the simulation was run forward in time to determine if Dance of the 
Planets achieved the simulation accuracy required to show the deflection of the asteroid 
and determine if the calculated deflection is sufficient to avoid collision. The results are 
presented in the next section. 

All MATLAB scripts, presented in the Appendix, were written in modular fashion. 
Running the script from the main program accesses additional script files as required. 

Each file can also be run from the Command window by typing the function name with the 
appropriate arguments provided. Help files are provided for each function file. The 
MATLAB scripts also plot the Earth and asteroid orbits. Aphelion, perihelion, the 
position of each body in their orbits is indicated on the plots. 
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7. Results and Discussion. 


Following the development of the MATLAB function files and orbital 
development routines, a problem was discovered with the accuracy of files imported to 
Dance of the Planets. Data files used to import elements into the simulation are space 
delimited and of the form : 

Name x y z Vx Vy Vz epoch 
The position coordinates were in units of E6 km and contained up to 5 decimal place 
precision. Velocity vector components were in km/sec and contained up to 6 decimal 
place precision. Epoch date was in Julian days and allowed inputs up to four decimal 
places. 

Importing a data file with this level of exactness would have been sufficient to 
accurately test the deflection hypothesis presented in section III. However, when data 
was imported, the values obtained were slightly different from the original file. Table 8 
show the result of importing the data file for asteroid 2340 Hathor. The file contained the 
consisted of the position and velocity vectors shown under Hathor, below, and was 
obtained by exporting actual position and velocity data from Dance of the Planets. The 
example column listed next to Hathor shows how the values differed following 
importation. 
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Ast. 

Hathor 

Example 

X 

13.29207 

13.29217 

y 

-75.41408 

-75.41404 

Z 

0.02352 

0.02351 

Vx 

44.735765 

44.735756 

Vy 

19.744673 

19.744721 

Vz 

-4.868622 

-4.868622 


Table 8. Imported Asteroid’s Altered Position and Velocity Vectors. 


The magnitude of the change between the desired values and the values as altered 
by Dance of the Planets is shown in Table 9. If the values had been imported correctly 
the position and velocity vectors would be identical. 


Asteroid 

Delta r (km) 

Delta v (m/s) 

Example 

108.1665 

0.1082 


Table 9. Position and Velocity Change from Original Orbit. 

While these errors were relatively small, they were of sufficient magnitude to 
impart significant errors at the initiation of the simulation and prevented precise and 
accurate testing of the theories; the magnitude of the change in velocity to be imparted to 
the asteroid was of approximately the same magnitude as the error generated by the 
program. Additionally, these errors were not consistent in magnitude or direction from 
one test to the next. The errors could not be biased out by altering the date of epoch for 
the simulation. 

Discussions with Terry Hannon of Applied Research & Consulting, Inc., one of 
the software developers for the program, lead to speculation that the problem may lie in 
the precision of the epoch currency. The program presents accuracies to one ten- 
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thousandth of a day, approximately eight seconds, while the accuracies required for the 
thesis testing was on the order of one second accuracy. In fact, although the epoch date is 
displayed out to four decimal place, the manufacturer believed the program accuracy to be 
somewhat less than this. 

The program is written in Pascal, and should the manufacturer make it available, 
the code could perhaps be modified so as to make it acceptable to conduct simulations of 
this type. Several simulations were conducted which did not require the importation of 
data files. These runs were very accurate and results were repeatable during multiple 
simulations. Overall, the program provides very accurate orbital positional data for 
predicting and visualizing solar system events. Without modification, however, the 
program is not sufficiently precise to conduct deflection simulations with changes of 
velocity on the order of 0.1 m/s. 

C. DEFLECTION RATIO SENSITIVITY TO ORBITAL PARAMETERS 

1. Introduction 

Since the original simulation program discussed in section B, above, was not 
satisfactory to test the deflection hypothesis, a simulation program was written in 
MATLAB which permitted the application of a given change in velocity to a variety of 
hypothetical asteroid orbits. The MATLAB scripts used are contained in Appendix C. 
The change in velocity was applied in the plane of the asteroid’s orbit. The direction of 
the change in velocity was applied in five equal increments from along the flight path 
direction to the anti-flight path direction. As the answers were symmetrically similar, the 
directional range was altered to five equally spaced inputs from along the flight path to 
orthogonal to the flight path. These are depicted in Figure 17. 
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Flight Path 



Figure 17. Deflection Directions With Respect to the Flight Path Vector. 

2. Variation in Orbital Period for a Circular Orbit 

The first sequence of tests applied a 0.135 m/s change in velocity to an asteroid 
with a circular orbit at 1 AU, inclined to the ecliptic. The Earth’s orbit was modeled as a 
circular orbit in the ecliptic at 1 AU. The simulation calculated the orbital elements of 
the asteroid’s orbits including the true anomaly. From those elements, the eccentric 
anomaly, mean anomaly, and mean motion were calculated. The orbital velocities at one, 
two, three, four, five, and ten years prior to collision were calculated, the desired change 
in velocity applied, and new orbital elements calculated. The entire procedure was then 
reversed. Finally, the asteroid’s position at the original time was determined, compared 
to the Earth’s position, and divided by the Earth’s radius to determine a miss ratio in 
Earth radii. As expected, the asteroid missed by two Earth radii for the flight-path and 
anti-flight path directions. Results show that the miss ratio decreased by the sine of the 
angle between the flight path and deflection direction. At 45 deg., the miss ratio 
decreased to 1.4 Earth radii, while at 90 deg., there was no deviation detected. The 
magnitude of the deflection distances were symmetric about 90 deg. 
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3. Constant Radius of Perihelion 


To demonstrate the accuracy of the mathematical solution developed in Chapter 
IV, miss ratios were calculated for an asteroid orbit with an fixed radius of perihelion of 1 
AU. Fixing perihelion distance results in varying eccentricity and semi-major axis 
lengths. Change in miss ratio with change in eccentricity and change in semi-major axis 
are shown in Figure 18 and Figure 19, respectively. Increasing eccentricity dictates an 
increase in semi-major axis for a constant radius of perihelion. As such, the two figures 
look very much the same. As the eccentricity increased from a minimum value of 0.0 to 
approximately 0.2, the miss ratio increased from 2 Earth radii to approximately 4.2 Earth 
radii. Over this same variation, the semi-major axis increased from 1 AU to 
approximately 1.25 AU. As shown in Eqn. (24), the deflection distance varies linearly 
with the asteroids period. But, because the asteroids period varies with a 312 , there is a 
relatively large increase in miss ratio for increasing values of eccentricity. 
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” Angle 

Figure 19. Miss Distance in Earth Radii vs. Deflection Angle vs. Semi-Major Axis. 

4. Variation in Eccentricity 

The effect of varying the eccentricity of the orbit while maintaining a constant 
semi-major axis length, shown in Figure 20, provided several interesting results. As the 
eccentricity increases and a is held constant, the perihelion distance decreases. The line 
of apsides must then rotate, increasing the true anomaly, to keep the collision distance at 
1 AU. The true anomaly increases approximately 5.7 deg. for each 0.1 increase in 
eccentricity. As the eccentricity increases, the miss distance decreases. The primary 
cause of this is the application of the change of velocity away from perihelion. The 
greatest deflections are achieved by applying the change in velocity at perihelion. The 
least efficient deflections are at aphelion. In this case, the change in velocity is applied 
approximately 91 deg. past perihelion for eccentricity of 0.01. This increases to 
approximately 102 deg. past perihelion at e = 0.2. The gradient of miss ratio with 
increase in eccentricity increases significantly beyond e = 0.2. This shows the importance 
of applying a given change in velocity at perihelion. 
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Figure 20. Miss ratio vs. Deflection Angle vs. Semi-Major Axes. 

5. Conclusions 

The results presented above are consistent with the mathematical conclusions 
presented in Chapter IV. One of the most important factors in maximizing the miss ratio 
is targeting the change in velocity to occur at perihelion. This is especially important for 
high eccentricity orbit with a semi-major axis close to that of the Earth. Much larger 
changes in velocity are required in these scenarios the further from perihelion the change 
is to occur. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. SUMMARY 

The issues surrounding ECA mitigation necessitate a wide variety of 
professionals, each with differing areas of expertise, form a project team to formulate the 
most feasible solution. This paper, for the Naval Postgraduate School and the Space 
Warfare Center, is the genesis of another small portion of that project team. Although in 
its infancy, this program, because of the unique composition of its members, can 
contribute unique and meaningful contributions to that effort. 

This paper serves several important purposes. First, it summarizes a few of the 
important issues involved in planetary defense and serves as a primer to helps to gain an 
understanding of ECAs. Secondly, it serves to develop and analyze the sensitivity of 
deflection angle and various orbital parameters to miss ratio. Most importantly, it 
identifies areas for further research and the continuation of this important problem. 

Dance of the Planets is not sufficient to serve as a model for the purpose of testing 
deflection hypotheses. The errors induced into orbital elements or vectors when 
importing data are of sufficient magnitude to prevent accurate analysis of results. Should 
the source code become available, increasing the precision used for calculations may 
solve this problem. This is not to say, however, that the program is without merit. Used 
as a visualization tool, it provides an extremely flexible and accurate model with which to 
observe intercept geometry. 

The variation of miss ratio with deflection angle, eccentricity, and semi-major axis 
changes was analyzed. Deflection along the flight path achieve the greatest miss ratios. 
Position in the orbit also play an important role in the magnitude of the miss ratio. For 
eccentric orbits, deflection at perihelion is important as miss ratio is fairly sensitive to 
true anomaly. As true anomaly at deflection increases, the miss ratio decreases. This a 






direct result of the deflection being applied at further angular distances from perihelion 
resulting in less of a change in period for each asteroid orbit. 

B. FURTHER RESEARCH OPPORTUNITIES 

There are three main areas in the asteroid mitigation problem which would benefit 
from additional research. First, a complete engineering and mathematical problem 
formulation should be accomplished. Expanding the scope of the current investigation to 
include a review of the state-of-the-art in change in velocity delivery options, intercept 
trajectory optimization, and weapons technology. 

Second, asteroid orbit determination accuracy versus the change in velocity 
required to mitigate a collision should be researched. The accuracy of ephemeris data is 
extremely important when calculating required changes in velocity and deflection angle. 
The relationship between the two should be determined to ensure that the proper 
mitigation steps are taken. 

Lastly, the formulation of the sensitivity model developed in this paper should be 
expanded to include gravitational perturbations and solve the problem using numerical 
integration techniques. A high fidelity model may provide additional insights into the 
sensitivity of miss ratio to deflection angle. 
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APPENDIX. MATLAB SCRIPT FILES 


This appendix contains the MATLAB function files constructed during the 
research of the asteroid mitigation problem. Section I contains the scripts which calculate 
a hypothetical asteroid’s position and velocity vectors based on the current position and 
velocity of the Earth. Demo.m calculates all data by accessing individual script files. In 
Demo.m the user inputs the Earth’s position and velocity vectors, rl and vl, respectively. 
Position is in 1E6 km and velocity is in km/sec. For an elliptical orbit, the user inputs into 
elip.m the desired semimajor axis length, eccentricity, and inclinatioa For a parabolic 
orbit, the user inputs into parab.m the desired semi-latus rectum, and inclination. 
Hyperbolic orbit script files were not developed. When demo.m is executed, available 
outputs include the asteroid’s position and velocity vectors in heliocentric ecliptic 
coordinates and the asteroid’s orbital elements. Current position in the orbit is output as 
true anomaly. 

Section 2 calculates the miss ratio, the ratio of the difference between the Earth’s 
position and the asteroid’s position one asteroid orbit following the imparted change in 
velocity to the radius of the Earth. The user inputs the Earth’s position and velocity 
vectors, as discussed above, into Demol.m. The user also inputs into demo.m the 
magnitude of the change in velocity to be applied, normalized for one orbit, and the 
number of orbits prior to collision the change in velocity is to be applied. In elipl .m the 
user inputs are the desired semi-major axis, eccentricity, inclination for asteroid at the time 
of Earth impact. The program then calculates the remaining orbital parameters for the 
asteroid’s orbit. From these, the eccentric anomaly, mean anomaly, and mean motion are 
found and used to find the asteroid’s position at some specified number of asteroid orbits 
prior to collision. The change in velocity is then applied in five 22.5 deg. increments from 
along the flight path to 90 deg. to the flight path. New values for mean motion, mean 
anomaly, and eccentric anomaly are found and the new position obtained at the original 
time. The miss ratio is then calculated to determine how many Earth radii the new 
position differs from the old position. Outputs include new velocity vectors, orbital 
elements, true anomaly, and miss ratio for each of the five directions. 
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Summary of function files: 

trunanom: 

True anomaly determination from r,v,e,n. 

incl: 

Inclination determination from h. 

lan: 

Long, of Ascending Node determination from n. 

aop: 

Argument of Perihelion determination from n,e,r. 

rot: 

Rotation matrix from from perifocal to HCI. Inputs are 

Omega,i, omega. 

elip: 

Finds r,v,nu,dcm from Earth’s rmag,dcm,nu. 

parabola: 

Finds r,v,nu,dcm from Earth’s rmag,dcm,nu. 

earthorb: 

Calculates Earth’s r vector in perifocal and HCI from a,ecc,dcm. 

cx,cy,cz: 

Compute a single axis direction cosine matrix. 

tlop: 

Compute true longitude of perihelion from e. 

deltav: 

Computes the new velocity vector components after the delta v is 
imparted to the asteroid. 

elements: 

Returns the matrix sets of: 

el__prev: orbital elements prior to deflection. 

el_next: orbital elements after deflection. 

rv_next: position and velocity vectors at original impact 

time. 

elipl: 

Computes the asteroid’s orbit given the Earth’s rmag,dcm,nu and 
user defined values for a,e,i. 

elip2: 

Computes the asteroid’s new r,v given the asteroid’s position at 
deflection, r, the dcm from the orbital elements, the position in the 
old orbit, ecc_a,a_a. 

functionl/2: 

Solve Kepler’s equation for E. 

param: 

Given r,v returns a,ecc,Omega,incl,omega,nu,E,M,dcm,rmag. 

Variables: 


Position vector. 

v: 

Velocity vector. 

e: 

Eccentricity vector. 

n: 

Node vector. 

h: 

Angular momentum vector. 

Omega: 

Longitude of the ascending node. 

i: 

Inclination. 

omega: 

Argument of perihelion. 

dcm: 

Direction cosine matrix from perifocal to heliocentric-ecliptic 
coordinate system for Earth’s orbit. 

nu: 

True anomaly. 

a: 

Semi-major axis. 

ecc: 

Eccentricity. 

dcm2: 

Direction cosine matrix from perifocal to heliocentric-ecliptic 
coordinate system for asteroid’s orbit. 
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*_old 

*_prev 

*_new 

*_next 

P 

r _P 

r_a 

E 

M 


Earth variable. 

Asteroid variable. 

Asteroid variable at time of impact. 

Asteroid variable at some time prior to impact. 

Asteroid variable just after application of change in velocity. 
Asteroid variable in new orbit but at time of original impact. 
Semi-latus rectum. 

Perihelion distance. 

Aphelion distance. 

Eccentric anomaly. 

Mean anomaly. 

Mean motion. 
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I. ORBIT DETERMINATION AND GRAPHING SCRIPTS 


A. CALCULATE AND PLOT ORBITS: FILE ENTITLED DEMO.M 


% Given: 

% First: 

% Second: 
% Third: 

% a. 

% b. 

% c. 

% Fourth: 
% a. 

% b. 


The Earth's position and velocity vectors. 

Calculate the Earth's orbital elements. 

Plot the Earth's orbit wrt the mean ecliptic. 

Choose the type of assailant object orbit: 

Circular or elliptical 

Parabolic 

Hyperbolic 

Choose the desired values for: 

Semi-major axis or parameter 
Eccentricity 

New inclination with respect to Earth's inclination. 


% Constants 


rtd=180/pi; 

dtr=l/rtd; 

I = [100]; 

J = [0 10]; 

K = [0 0 1]; 

AO = 1.4959965e8; 
TO = 5.0226757e6; 
SO = 29.784852; 


% Radians to degrees 
% Degrees to radians 
% Onit vectors 


% km 
% sec 
% km/sec 
% km A 3/sec*2 


% Earth values: 

rl= [-26 144.8 -.00038]; % Position 

vl= [-29.8 -5.376 -.000043]; % Velocity 


r=rl*le6/AO; 

v=vl/SO; 

rmag=norm(r); 
vmag=norm(v); 


hmag=norm(h) ; 

% eccentricity vector 

e=cross (v,h) -r/rmag,- 

if ecc<le-10, % Check if eccentricity is zero 
e = [0 0 0]; 
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% node vector 

nmag=norm(n) ; 

if nmag<le-10, % Check if inclination is zero 
n = [0 0 0]; 
nmag=0; 


% Find the true anomaly at epoch 
nu = truanom(r,v,e,n) ; 

% Compute the inclination 
inc = incl(h) 

% Compute the longitude of ascending node 
Omega=lan(n) 

% Compute the argument of perigee 
omega = aop(n,e,r) 

% Compute the rotation matrix 
dcm=rot(Omega,inc,omega) ; 

V Determine if user wants another orbit calculated 
pit = input('Do you want to calculate a new orbit? y/n [y]:',’s '); 

if isempty(plt) 
pit = 'y'; 


if pit == 'y', 


% choose an a,e for orbit and solve for nu_2 for the new orbit 

fprintf( 1 \nYou may now calculate a new intercept orbit type.\n'); 

type = input('Do you want to plot an Ellipse, Parabola, or Hyperbola? e/p/h 


[e] : 


■ '); 


if type == 'e', % Eliptical/Circular Plot 

[r_new,v_new,nu_new,dcm2]=elip(rmag,dcm,nu); 

elseif type == 'p', % Parabolic Orbit Plot 

[r_new,v_new,nu_new]=parabola(rmag,dcm,nu) 

elseif type == 'h', % Hyperbolic Orbit Plot 

[r_new,v_new,nu_new]=hyperb(rmag,dcm,nu) 
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hold 

plot: 


:,yy,zz. 


xyz_r=r; % position 

plot3(xyz_r(l),xyz_r(2) ( xyz_r(3),'g+'); 

if type == 'e' , % Eliptical/Circular Plot 

legend(' b' , ' Asteroid' 1 s 

orbit' , 'or','Perihelion','*r','Aphelion','+g','Current Position',... 
'm','Earth''s orbit') 

elseif type == 'p' | 'h', % Parabolic Orbit Plot 

legendt'b','Asteroid''s orbit','or','Perihelion',’+g','Current 
Position'.... 

'm','Earth''s orbit') 


xyzj>er=dcm* [rp;0;0] ; % perihelion 

plot3 (xyz_per (1) ,xyz_per (2) ,xyz_per(3) , 'mo') ; 

xyz_apo=dcm*[-ra;0;0]; % aphelion 

plot3(xyz_apo(1),xyz_apo(2),xyz_apo(3),'m*'); 

% end plotting of' second orbit 
% Axis plotting statements. 

1=2; 
ds =.2; 

A= [ 0 0 0 ; 1 0 0] ; 

B= [0 0 0; 0 1 0] ; 

C=[0 0 0; 0 0 1] ; 
plOt3(A,B,C,'w'); 

text (1+ds, 0,0, 'X' , 'horizontalalignment' , ' center') ,- 
text(0,1+ds,0,'Y','horizontalalignment','center'); 
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text (0,0,1+ds,'Z', 1 horizontalalignment','center'); 

view(viewmtx(135,30,25)); 

title('General View'); 

axis('square');axis off; 

axis ( [-a,a,-a,a,-a, a]) ; 


% View radio buttons 
txt_view=uicontrol(gcf,... 

'Style','text',... 

'String','Asteroid Orbit View',... 

'Units','normalized' , . . . 

'Position', [.05,.87,.20,.03] ) ; 
view_gen=uicontrol(gcf,... 

'Style','radio',... 

'String','Standard',.. . 

'Units',’normalized',. .. 

'Position', [.05, .80, .20,.05]_ 

'Value',1,... 

1 CallBack',[... 

' title(''General View''),'... 

' set(view_gen,''Value'',1),'... 

' set(view_x,''Value'',0),'... 

’ set(view_y,''Value'',0),'... 

'set(view_z,''Value'',0),'... 

' set(view_n,''Value'',0),'... 

'view (135,30) '] ) ; 
view_x=uicontrol(gcf, . . . 

'Style','radio', . .. 

'String','Down X Axis', ... 

'Units','normalized'.... 

'Position', [.05,.75,.20,.05] , . . . 

'CallBack',[... 

'title(''View Looking Down +X axis''),'... 

'set(view_gen,''Value' ',0) , ' . . . 

'set(view_x,''Value'',1) , ' ... 

'set(view_y,''Value'',0) , ' ... 

'set(view_z,''Value'',0) , ' . .. 

'set(view_n,''Value'',0), ' ... 

•view( [1,0,0])']); 

view_y=uicontrol(gcf, . . . 

'Style','radio' .... 

'String','Down Y Axis',... 

'Units','normalized' , . . . 

'Position', [.05,.70,.20,.05]_ 

'CallBack',[... 

'title(''View Looking Down +Y axis''),'... 
'set(view_gen,''Value'',0) , ' ... 

’set(view_x,''Value'',0) , ' ... 

'set(view_y,''Value' ' , 1) , ' . . . 

'set(view_z,' 'Value'',0) , ' . . . 

'set(view_n,''Value' ' , 0) , ’ . . . 

'view( [0,1,0] )']) ; 
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view_z=uicontroi(get, . . . 

'Style','radio'.... 

'String','Down Z Axis' , ... 

'Units','normalized'.... 

'Position', [.05, .65,.20,.05]_ 

'CallBack',[... 

'title(''View Looking Down +Z axis''),'... 

'set(view_gen,''Value'' , 0) , ' . . . 

'set <view_x,''Value'',0),'... 

'set(view_y,''Value'',0),'... 

'set(view_z,''Value'',1),'... 

'set(view_n,''Value'',0),'... 

'view( [0,0,1]) ']); 

xyz_norm=dcm2 * [ 0; 0; 1 ] ,- 
view_n=uicontrol(gef, ... 

'Style','radio' , .. . 

'String','Down Orbit Normal',... 

'Units','normalized',... 

'Position', [.05, .60,.20,.05],... 

’CallBack’,[... 

'title( 1 ’View Looking Down Orbit Normal•'),'... 

'set(view_gen,''Value'' , 0),'... 

'set(view_x,''Value'' , 0) , ' ... 

'set(view_y,''Value'',0),'... 

1 set(view_z,’’Value 11 ,0),’... 

1 set(view_n, 1 ’Value'',1),'... 

'view( [xyz_norm(1) ,xyz_norm(2) ,xyz_norm(3) ])']),- 


% Zoom Slider 
zoom=uicontrol(gef, . . . 

'Style','slider' , . . . 

'Units','normalized',... 

'Position',[.05,.15,.20,.05],... 

'Min',.5,'Max',5.5,'Value',1_ 

'CallBack',[... 

'set(zoom_cur,''String'',num2str(get(zoom,''Val 
' set (gca,' 'Xlim' ’,1- 

get (zoom,''Value'’),get(zoom,’’Value 1 ’)],’,... 

’ 1 ’Ylim'',[-get(zoom,''Value''),get(zoom,''Value'')],',... 
'''Zlim'',[-get(zoom,''Value''),get(zoom,''Value'')])']); 
zoom_label=uicontrol(gef,... 

'Style','text',... 

'Units','normalized',... 

'Position', [.05, .205,.20, .03).... 

'String','Zoom'); 
zoom_cur=uicontrol(gef,... 

'Style','text',... 

'Units','normalized’ , . . . 

'Position',[.095,.115,.11,.03] , ... 

' String' ,num2str (get (zoom, 'Value') ) ) ,- 
zoom_in=uicontrol(gef, . .. 

’Style','text' , ... 

'Units','normalized' , . . . 


75 















'Position', [.05,.115, 
* String’,'In’); 
zoom_out=uicontrol(gcf, . . . 
'Style','text' , . . . 

1 Units','normalized' , 




button 


B. AOP 



Calculation of the argument of 
perigee, in degrees 

LCDR Wade Knudson 

Naval Postgraduate School 

Orbit Visualization Routines 

Inputs: n,e,r 

Output: omega 

Files called: tlop 

[omega] = aop(n,e,r) 


C. ARGUMENT OF PERIHELION 


mag=norm(n); 

mag=norm(r); 
= [10 0]; 

= [ 010 ]; 

: = [0 0 1 ]; 
td=180/pi; 






% Check to see if the orbit is elliptical (e>0) and in the ecliptic (n=0) 
% if so, then omega is not used; use true longitude of periapsis (Vallado 
p. 136) 

elseif nmag == 0 & ecc ~=0, 
om_true = tlop(e); 

%fprintf('This is an elliptical orbit in the ecliptic.\n') 
%fprintf('The inclination is zero.\n') 

%fprintf('\n') 

Srfprintf ('The angle between Aries and periapsis is called \n') 
%fprintf('the true longitude of periapsis, and is %3.1f 


% Check to see if the orbit is a circular (e=0) and inclined to the 
ecliptic 

elseif nmag ~= 0 & ecc == 0, 
if dot (r,K) <0, 

omega = 360 - acos(dot(n,r)/(nmag*rmag))*rtd; 




omega = acos(dot(n,r)/(nmag*rmag))*rtd; 


% All remaining cases 
else 

if dot (e,K) < 0, 

omega = 360 - acos(dot(n,e)/(nmag*ecc))*rtd; 
else 

omega = acos(dot(n,e)/(nmag*ecc))*rtd; 


D. CX 

function f=cx(angle) 

f=[l 0 0; 0 cos(angle) sin(angle); 0 -sin(angle) cos (angle) ] ,- 

E. CY 

function f=cy(angle) 

f= [cos (angle) 0 -sin(angle) ,- 0 1 0; sin(angle) 0 cos (angle) ] ; 

F. CZ 

function f=cz(angle) 

f=[cos(angle) sin(angle) 0; -sin(angle) cos(angle) 0; 0 0 1]; 
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G. EARTH ORB 


% 

% Calculation of the Earth's orbit 


% 

% LCDR Wade Knudson 

% Naval Postgraduate School 

% Orbit Visualization Routines 

% 

% Inputs: a,ecc,dcm 

% Output: x,y,z: perifocal coordinates 

% xx,yy, zz: heliocentric coordinates 

% 

% Files called: None 

% 

% [x,y,z,xx,yy,zz] = earthorb(a,ecc,dcm) 

% 

% --- 


function [x,y,z,xx,yy,zz] = earthorb(a,ecc,dcm) 


m=400; 

E=linspace(0,2*pi,m) 

x=-a*ecc+a*cos (E) ;y=a*sqrt (l-ecc A 2) *sin(E) %Battin p. 158 

z=zeros(l,m); 
for k=l:tn 

point=dcm* [x(k) ;y(k) ;z(k)] ; 

xx (k) =point (1) ;yy (k) =point (2) ;zz (k) =point (3) ; 

end; 


H. ELIPTICAL ORBIT CALCULATIONS 


Calculation of the Eliptical/circular orbit 


LCDR Wade Knudson 

Naval Postgraduate School 

Orbit Visualization Routines 

Input s: rmag,dcm,nu 
User inputs: a,ecc,incl 
Output: r,v, and nu for the new orbit 

Files called: cx,cz 
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of node 




inclination, 

% and finally rotate nu_2 to get epoch to r. 

% Asteroid's dcm matrix 

dcm2= dcm*cz(-nu*dtr)*cx(-newinc*dtr)*cz(nu_new*dtr) 


m=200; 

E=linspace(0,2*pi,m) ; 
x=-a*ecc+a*cos(E) 

y=a*sqrt(l-ecc A 2)*sin(E); % Battin p. 158 

z=zeros (1 ,tn) 
for k=l:m 

point=dcm* [x (k) ;y (k) ; z (k) ] ; 

xx (k) =point (1) ;yy(k) =point (2) ,-zz (k) =point (3) ; 
point2=dcm2*[x(k);y(k);z (k)]; 

xxx(k) =point2 (1) ;yyy (k) =point2 (2) ,-zzz (k) =point2 (3) ; 
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% Plot the asteroid's new orbit 
hold on; 

plot3(xxx,yyy,zzz, ’b’) ; % orbit 

xyz_per=dcm2*[rp;0;0] ; % perihelion 

plot3(xyz_per(l),xyzjper(2),xyz_per(3),'bo'); 

xyz_apo=dcm2*[-ra;0;0]; % aphelion 

plot3(xyz_apo(l),xyz_apo(2),xyz_apo(3),'b*'); 

xyz_r=(dcm2*r 1 )'; % position 

plot3(xyz_r(1),xyz_r(2),xyz_r(3) ,'g+'); 
legend('Asteroid' 's orbit','Perihelion','Aphelion') 

% Calculate the new r,v vectors 
r = xyz_r; 
v = (dcm2*V ) » ; 


I. INCLINATION 

%---- 

% 

% Calculation of the inclination in degrees 

% 

% 

% LCDR Wade Knudson 

% Naval Postgraduate School 

% Orbit Visualization Routines 

% 

% Inputs: h 

% 

% Output: inc 

% 

% Files called: none 

% 

% [inc] = inc (h) 


function inc = inc(h) 

K =[0 0 1 ] ; 
rtd=180/pi; 
hmag=norm(h); 


inc 


(dot(h,K)/hmag)*rtd; 






J. LONGITUDE OF ASCENDING NODE 


.se % If some inclination t 

% find Omega 

if dot (n, J) < 0, 

Omega = 360 - acos(dot(n,I)/nmag)*rtd; 


LCDR Wade Knudson 

Naval Postgraduate School 

Orbit Visualization Routine 


I 









% Files called: cx,cz 

% [r,v,nu_new] = parabola(rmag,dcm,nu) 

% 

% -- 

function [r,v,nu_new] = parabola(rmag,dcm,nu) 

% Constants 

rtd = 180/pi; 
dtr = 1/rtd; 


nu_new=acos((p/rmag-1)/ecc)*rtd; % true anomaly 

rp = p/(l+ecc) % perihelion 

% Perifocal Coordinates 

r = [rmag*cos(nu_new*dtr) rmag*sin(nu_new*dtr) 0]; 
v = sqrt(mu/p)*[-sin(nu_new*dtr) ecc+cos(nu_new*dtr) 03; 


% rotate the new orbit to the line of node (-nu), then to desired 
inclination, 

% and finally rotate nu_2 to get epoch to r. 

dcm2= dcm*cz(-nu*dtr)*cx(-newinc*dtr)*cz(nu_new*dtr); 


m=200; 

f-linspace(-pi+.l,pi-.l,m) ; 
xl » (p*cos (f)) ;yl = (p*sin(f)); 
x2 = (1+cos (f));y2 = (l+cos(f)); 

x = xl./x2,-y = yl./y2; % Battin p. 158 

z=zeros(l,m); 

for k=l:m 

point=dcm* [x(k) ;y(k) ;z(k) ] ; 

xx (k) =point (1) ;yy (k) =point (2) ;zz (k) =point (3) 
point2=dcm2* [x(k) ; y(k) ; z (k) ] ; 

xxx(k) =point2 (1) ;yyy(k) =point2 (2) ,-zzz (k) =point2 (3) ; 


% Plot the asteroid's new orbit 
hold on; 
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p!ot3(xxx,yyy,zzz,'b'); % orbit 

xyz_per=dcm2* [rp,-0;0] ; % perihelion 

plot3(xyz_per(1),xyz_per(2),xyz_per(3),'bo'),- 

xyz_r=(dcm2*r’)'; % position 

% rvector 

plot3(xyz_r{1),xyz_r(2),xyz_r(3),'g+'); 

% Calculate the new r,v vectors and orbital elements, 
r = xyz_r 
v = (dcra2*v') ' 


L. PRINT 

% Puts a print push button on a figure, 
function prnbut = prt 

prt=uicontrol(gcf, .. . 

'Style','push',... 

'Units','normalized' , . . . 
'Position', [.05,.02,.20, .05], . . . 
1 String 1 , 1 Print',... 

'CallBack','print'); 


M. ROTATION DCM 


% 

% Calculation of the rotation matrices for a 3-1-3 
% rotation to convert from perifocal coordinates to 
% heliocentric-ecliptic 

% 

% LCDR Wade Knudson 

It Naval Postgraduate School 

% Orbit Visualization Routines 

% Inputs: Omega,inc,omega 

% 

% Output: dcm 

% 

% Files called: cx,cz 


% [dcm] = rotmat(Omega,inc,omega) 

%--- 


function [dcm] = rot(Omega,inc,omega) 


dtr=pi/180; 

0=0mega*dtr ,- 
i=inc*dtr ; 
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o=omega*dtr; 

dcm = cz (-0) *cx(-i) *cz (-o) ; 

N. TRUE LONGITUDE OF PERIHELION 

%- 

°s 

% True longitude of periapsis 

% Used for elliptical orbits in the reference plane 

% LCDR Wade Knudson 

% Naval Postgraduate School 

% Orbit Visualization Routines 

V 

% Output: True longitude of periapsis 
% 

% Files called: none 

% 

% [om_true] = tlop(e) 

% 

% ------- 


function om_true = tlop(e) 

rtd = 180/pi; 

I = [10 0]; 

J = [0 10]; 
emag = norm(e); 

if dot (J, e) < 0, 

om_true = 360 - acos(dot(I,e)/emag)*rtd; 

else 

om_true = acos(dot(I,e)/emag)*rtd; 

O. TRUE ANOMALY 


% 

% Calculation of the true anomaly at epoch, in degrees 
% 

% LCDR Wade Knudson 

% Naval Postgraduate School 

%• Orbit Visualization Routines 

% 

% Inputs: r,v,e 

% Output: nu 

% 


Files called: 






function [nu] = truanom(r,v,e,n) 


nmag=norm(n); 

rmag=norm(r); 

I = [10 0]; 

J = [0 10]; 

K = [0 0 1]; 
rtd=180/pi; 

% First check for circular orbit in the ecliptic. If not, skip, 
if ecc == 0 & n == 0, 
if dot(J,r)<0, 

lamda_t = 360 - acos(dot(I,r)/rmag)*rtd; 
else 

lamda_t = acos(dot(I,r)/rmag)*rtd; 
nu = lamda_t; 

%fprintf('This is a circular orbit in the ecliptic. \n') 

%fprintf('The longitude of ascending node is undefined.\n') 

%fprintf('The inclination is zero.\n') 
fprintf('\n') 

%fprintf('The angle between aries and the position of the body\n') 
%fprintf('is called the true longitude at epoch, and is %3.lf deg.\n',nu) 


% Next check for circular inclined orbit. If not, skip, 
elseif ecc == 0 & nmag -= 0, 

if dot (r,K) <0, 

u = 360 - acos(dot(n,r)/(nmag*rmag))*rtd; 
else 

u = acos (dot (n, r) / (nmag*rmag)) *rtd,- 


fprintf('This is a circular, inclined orbit. \n') 
fprintf('\n') 

fprintf('The argument of perigee is known as the \n') 
fprintf ('argument of latitude at epoch, and is %3.1f deg.\n',u) 

else 

if dot(r,v) < 0, 

[nu] = 360 - acos(dot(e,r)/(ecc*rmag))*rtd; 
else 

[nu] = acos(dot(e,r)/(ecc*rmag))*rtd; 
end 

%fprintf(* This is an elliptical orbit. \n') 
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if nmag == 0 & ecc ~= 0, 

%fprintf('The orbit is in the ecliptic.\n') 

else 

%fprintf('The orbit is inclined to the eclipt 
%fprintf('\n' ) 

%fprintf('The true anomaly at epoch is %3.1f deg.\n',nu) 




II. SENSITIVITY SCRIPTS 


A. MISS RATIO CALCULATIONS: DEMOl 

clear all 
rtd=18 0/pi; 
dtr-l/rtd; 

I = [100]; 

J = [0 10]; 

K = [0 0 1]; 

AU = 1.4959965e8; % km 

TO = 5.0226757e6; % sec 

SO = 29.784852; % km/sec 

mu = 1; % km*3/sec A 2 

%r=[-26 144.8 -.00038]; 

%v=[-29.8 -5.376 -0.000043]; 


% Calculate the orbital elements for the Earth from the r/v vectors,- 
% Input are the position and velocity vectors 
% Outputs are shown: 

[a_e,ecc_e,0_e,inc_e,o_e,nu_e,E_e,M_e,dcm_e / rmag_e] = param(r.v) ; 

% 

% choose an a,e for orbit and solve for nu_2 for the new orbit 


% Determine if user wants another orbit calculated 
%plt = input('Do you want to calculate a new orbit? y/n [y]:','s'); 

%if isempty(plt) 
pit = 'y'; 

%end 

if pit «= 'y', 

%fprintf('\nYou may now calculate a new intercept orbit type.\n') ; 

%type = input('Do you want to plot an Ellipse, Parabola, or Hyperbola? 
e/p/h [e] :','s'); 

%if isempty(type) 
type = 'e'; 

%end 
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% Calculate new r/v vectors based on Inputting desired values of 
% a,e,inc. Omega,omega, and nu are not variables, but fixed, based 
% on the Earth's orbit. 

if type == 'e', % Eliptical/Circular Plot 

[r_new,v_new,nu_new,dcm2,ecc]=elipl(rmag_e,dcm_e,nu_e); 

elseif type == 'p', % Parabolic Orbit Plot 

[ r_new,v_new,nu_new]=parabola(rmag,dcm,nu) 

elseif type == 'h', % Hyperbolic Orbit Plot 

[ r_new,v_new,nu_new]=hyperb(rmag,dcm,nu) 

end % end plotting of second orbit 


% Compute the orbital parameters of the new orbit: 
global compare 

[a_a,ecc_a,0_a,inc_a,o_a,nu_a,E_a,M_a,dcm_a,rmag_a] = param(r_new,v_new) ; 

% Compute the period of the asteroid’s orbit 
mu_s = 1; 

TP = 2*pi*(a_a)"l.5/sqrt(mu_s)*TU/3600/24/365.25; % In Earth years 


% Input here how many asteroid orbits to go through. You should go through 
% at least one asteroid orbital period to ensure sufficient time for 
deflection. 

n_a= sqrt(mu_s/a_a A 3) ; % mean motion 

t=-TP*num_orbs*365.25*24*3600/TU; % convert years to canonical 

global t 

% Compute the required delta V required along the flight path 
del_tot=0.000135/num_orbs; 

% Compute mean anomaly for some past time. 

M_prev = M_a + n_a*t; 

M_prevl = M_prev - fix(M_prev/2/pi) *2*pi + 2*pi; 
global M_prev ecc_a 


% Previous eccentric anomaly 
E_prev = fzero('fun',M_prev) ; 

% Compute previous true anomaly 
nu_prev = acos((ecc_a-cos(E_prev))/(ecc_a*cos(E_prev)-1))*rtd; 

% Compute rmag for previous time 
rmag_prev = a_a*(l-ecc_a*cos(E_prev)) ; 

% Compute the r/v at the previous time 










[r_prev v_prev]=elip2(rmag_prev,dcm_a,nu_prev,ecc_a,a_a); 

% [ a_a_p, e c c_a_p, 0_a_p, i nc_a_p, o_a_p, nu_a_p, E_a_p, M_a_p, dcm_a_p ] = 
param(r_prev,v_prev) 


% Compute the new matrix of delta_v's 

tel_prev el_next rv_next]=deltav(r_prev,v_prev,del_tot,nu_new) ; 

delta_rv=zeros(6,6); 
for s=l:6 

delta_rv(s,1:3) = (r_e - rv_next(s,1:3))*le6; 
delta_rv(s,4:6) = v_prev - rv_next(s,4:6); 

mag_del_rv(s,:) = [norm(delta_rv(s, 1 :3)) norm(delta_rv(s,4:6))]; 

end 

time_ratio= [nu_new a_a ecc_a (mag_del_rv(2:6,1)/6378.145)']; 
fprintf ( 1 \n nu a e 0.0 22.5 45 67.5 

90\n 1 ); 
fprintf( 1 \n 

%7.4f%7.4f%7.4fSr7.4f%7.4f%7.4f%7.4f%7.4f%7.4f\n\n',[time_ratio]); 

fprintf(’\n') ; 
keyboard 

B. DELTA V DETERMINATION 

function [el_prev,el_next,rv_next] = deltav(rl,vl,dv_tot,nu_new) 

% Constants 
rtd=180/pi; 
dtr=l/rtd; 

I = [10 0]; 

J = [0 10]; 

K = [0 0 1]; 

AU = 1.4959965e8 
TU = 5.0226757e6 
SU = 29.784852; 


% Orbit values: 
r=rl*le6/AU; 
v=vl/SU; 

vmag=norm(v); 

% angular momentum 
h=cross(r,v); 

hmag=norm(h) ; 

% eccentricity vector 
e=cross(v,h)-r/rmag; 
ecc=norm(e); 

if ecc<le-10 ( % Check if eccentricity is zero 

e - [0 0 0]; 
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% km/sec 
% km' , 3/sec*2 









% node vector 
n = cross (K,h) ,- 

nmag=norm(n) ; 

if nmag<le-10, % Check if inclination 

n = [0 0 0],- 
nraag=0,- 

% Find the true anomaly at epoch 
nu = truanom(r,v,e,n); 

% Compute the inclination 
inc = incl (h) ,- 

% Compute the longitude of ascending node 
Omega=lan(n),- 

% Compute the argument of perigee 
omega = aop(n,e,r); 

% Compute the rotation matrix 
dcm=rot(Omega,inc,omega); 

direct=linspace(pi/2,0,5); 
x=cos(direct); y=sin(direct) ; 
z=zeros(1,5); 


for k = 1:5, 

dv= [x (k) y (k) z(k)]; 

if abs(x(k))<le-12 
x(k)=0; 

if abs(y(k))<le-12 
y(k)=0; 

end 

if abs(y(k))<le-12 
y(k)=0; 


dvl = (dv_tot* [x(l) y (1) 
dv2 = (dv_tot* (x(2) y(2) 
dv3 = (dv_tot*[x(3) y(3) 


z(l)] ') 
z(2) ] ') 
z (3) ] ') 
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dv4 


= (dv_tot* [x(4) y (4) z (4) ] ’) ; 
dv5 = (dv_tot* [x (5) y(5) z(5)]'); 

rot = cz (nu_new*dtr) ,- 

(rot*(vl + dvl 
(rot*(vl + dv2 
(rot*(vl + dv3 
(rot*(vl + dv4 
(rot*(vl + dv5 



% Print the old position and velocity vectors c 




r_prev = rl; 

dv = [v_prev;vnew_l,-vnew_2;vnew_3;vnew_4,-vnew_5] ; 


%fid=fopen('new_orb.vec','w'); 

%fprintf(fid,'ObjVectors x y z vx 

vz epoch\n'); 

%fprintf(fid,'%6.Of%ll.5f%11.5f%ll.5f%10.6f%10,6f%10.6fSrll.4f\n', 
v_prev epoch] ) ; 

%fprintf(fid, '%6.Of%ll.5f%ll.5f%ll.5f%10.6f%10.6f%10.6f%ll,4f\n' , 
vnew_l epoch]); 

%fprintf(fid,'%6.Of%ll.5f%ll.5f%ll.5f%10.6f%10.6f%10,6f%11.4f\n', 
vnew_2 epoch]); 

%fprintf(fid,'%6 . Of%ll.5f%ll.5f%ll,5f%10.6f%10.6f%10.6f%ll,4f\n', 
vnew_3 epoch]); 

%fprintf(fid, '%6.Of%ll.5f%ll.5f%ll.5f%10.6f%10.6f%10.6f%ll.4f\n', 
vnew_4 epoch]); 

%fprintf(fid, '%6.Of%ll.5f%ll.5f%ll.5f%10,6f%10.6f%10,6f%11,4f\n', 
vnew_5 epoch] ) ,- 


vy 

[0 r_prev 
[1 r_prev 
[2 r_prev 
[3 r_prev 
[4 r_prev 
[5 r_prev 


%fprintf(•\nObjVectors x y z vx 

%fprintf('%6.0f%11.5f%11.5f%11.5f%10.6f%10.6f%10.6f\n’,[0 
%fprintf('%6.0f%11.5f%11.5f%11.5f%10.6f%10.6f%10.6f\n',[1 
%fprintf(’%6.0f%11.5f%11.5f%11.5f%10.6f%10.6f%10.6f\n',[2 
%fprintf(’%6.0f%ll.5f%ll.5f%ll.5f%10.Sf%10.6f%10.6f\n',[3 
%fprintf(•%6.0f%11.5f%11.5f%11.5f%10.6f%10.6f%10.6f\n 1 , [4 
%fprintf(’ * 6.0 f£ll.5f%ll.5f%ll.5f%10.6f%10.6f%10.6f\n',[5 


vy 

r_prev v_prev ]) 
r_prev vnew_l ]) 
r_prev vnew_2 ]) 
r_prev vnew_3 ]) 

r_prev vnew_5 ]) 


(el_prev,el_next,rv_next]=elements(rl,dv,epoch,nu_new); 

C. DETERMINATION OF ORBITAL ELEMENTS 

% First bring in the r and v vectors that you want to convert. 

% r is the position vector 

% dv is a 5x3 matrix with the velocity vectors for the various delta v's 
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:] = elements (r,dv, epoch, 


. [el_prev,el_next,rv_ne 


% Constants 
rtd=180/pi; 
dr.r«l/rtd; 

I » [10 0] ; 

J ~ [0 1 0] ; 

K 1.0 0 1] ; 
at; -- j.. 4959965e8 ; 
TO - 5.0226757e6; 
SU -■ 29.784852; 


% km 
% sec 
% km/sec 
% km*3/sec*2 


% Orbit values; 

[m n] =size (dv) ; 
for s = l:m. 


% Convert to canonical units 
r=rl*le6/AU; 
v=vl(s, :)/SU; 

rmag=norm(r); 
vmag=norm(v) ; 

% angular momentum 

% eccentricity vector 
e=cross(v,h)-r/rmag,- 
ecc_n=norm(e); 

if ecc_n<le-10, % Check if eccentricity is zero 

e = [0 0 0]; 
ecc_n=0; 

% semi-major axis 
p = hmag^2/mu; 
a = p/(l-ecc_n^2); 

% node vector 
n = cross (K,h) 

nmag=norm(n); 

if nmag<le-10, % Check if inclination is zero 

n - [0 0 0]; 
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nmag=0; 

% Find the true anomaly at epoch 
nu = truanom(r,v,e,n) ,- 


% Compute the inclination 
inc = incl(h); 

% Compute the longitude of ascending node 
Omega=lan(n); 

% Compute the argument of perigee 
omega = aop(n,e,r); 

% Compute the rotation matrix 
dcm=rot (Omega, inc,omega) ,- 

% Compute the eccentric anomaly and mean anomaly 
E = acos ( (ecc_n + cos(nu*dtr)) / (1+- ecc_n*cos(nu*dtr)) ); 

if nu>180 

E=2*pi-E; 


M = (E - ecc_n*sin(E)) ; 

H=0; 

n = sqrt(mu_s/a*3); 
global t 

M_next = M - n*t; 

M_nextl = M_next - fix(M_next/2/pi)*2*pi; 
global M_next ecc_n 

% Next eccentric anomaly 
E_next = fzero('funl',M_next); 

% Compute next true anomaly 

nu_next =•acos((ecc_n-cos(E_next))/(ecc_n*cos(E_next)-1))*rtd; 

if sin(E_next)<0, 

nu_next =3 6 0-nu_next; 

end 

% Compute rmag for next time 
rmag_next = a* (l-ecc_n*cos (E_next) ) ,- 

% Compute the r/v at the next time 

[ r_next v_next ] =e 1 ip2 (rmag_next, nu_new, nu_next, ecc_n, a) ,- 



el_prev(s,:)=[a ecc_n Omega inc omega M n E]; 
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el_next(s,:) = [M_nextl E_next nu_next rmag_next] ; 
rv_next(s,:)=[r_next v_next]; 




D. ELIP 1 

% --- 

% 

% Calculation of the Eliptical/circular orbit 

% 

% LCDR Wade Knudson 

% Naval Postgraduate School 

% Orbit Visualization Routines 

% 

% Inputs: rmag,dcm,nu 

% User input: a,ecc,newinc 

% 

%• Output: r,v, and nu for the new orbit 
% 

% Files called: cx,cz 

% [r,v,nu_new] = elip (rmag,dcm,nu) 

% 

%-- - 


function [r,v,nu_new,dcm2,ecc] = elipl(rmag,dcm,nu) 


% Constants 
rtd = 180/pi; 
dtr = 1/rtd; 

AU = 1.4959965e8; 
TU = 5.0226757e6; 
SU = 29.784852; 


H km 
% sec 
% km/sec 
% km A 3/sec*2 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


% These are the controlling parameters of the new ellipse 
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


r_p=a*(1-ecc); 

p= r_p*(1-ecc); %AO /(1-ecc)/AO; 


% If circular then no 
% rotation to argument 
% of periapsis (nu_new) 
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ir penapsis 
s apoapsis 


%Perifocal coordinates 

rl = [rmag*cos(nu_new*dtr) rmag*sin(nu_new*dtr) 0]; 

vl = sqrt(mu/p)*[-sin(nu_new*dtr) ecc+cos(nu_new*dtr) 0]; 

h=cross(rl,vl); 
e=cross(vl,h)/mu-r/rmag; 


% rotate the new orbit to the line of node (-nu), then to desired 
% and finally rotate nu_2 to get epoch to r. 

dcm2= dcm*cz (-nu*dtr) *cx (-newinc*dtr) *cz (nu_new*dtr) 
xyz_r=(dcm2*rl') 1 ; 

r = xyz_r*AU/le6; 
v = (dcm2*vl')'*SU; 


E. ELIP 2 

%...-.. 

% 

% Calculation of the Eliptical/circular orbit 
% 

% 

% LCDR Wade Knudson 

% Naval Postgraduate School 

% Orbit Visualization Routines 

% 

% Inputs: rmag,dcm,nu 

% Output: r,v, and nu for the new orbit 

% 

% Files called: cx,cz 
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[r,v,nu_new] = elip(rmag,dcm,nu) 


% 

% - 

function [r,v] = elip2(rmag,dcm,nu_prev,ecc,a) 

% Constants 
rtd=180/pi; 
dtr-l/rtd; 

X = [10 0]; 

J = [010]; 

K = [0 0 1]; 

AU = 1.4959965e8; % km 

TU =. 5.0226757e6; % sec 

SO = 29.784852; % km/sec 

mu = 1; % km' , 3/sec > ‘2 

p = a*(l-ecc A 2); 
rp = p/(1+ecc); 
ra = p/(1-ecc); 

temp = cz(nu_prev*dtr); 

r = (temp*[rmag*cos(nu_prev*dtr); rmag*sin(nu_prev*dtr); 0])'; 
v = (temp* (sqrt (mu/p) ) * [-sin (nu_prev*dtr) ; ecc+cos (nu_prev*dtr) 0] ) ' ; 

r=r*AU/le6; 
v=v* SU; 


% periapsis 
% apoapsis 


F. FUNCTION 

% Solve Kepler's equation for E 
function x = fun(E) 
global M_prev ecc_a 

x = M_prev - E + ecc_a*sin(E) j 

G. FUNCTION 1 

% Solve Kepler's equation for E 
function [x] = funl(E) 
global M_next ecc_n 


x = M_next - E + ecc_n*sin(E) 

H. PA RAM 

% Solves for a,ecc,Omega,i,omega,nu,E,M,dcm, and rmag from r,v 

function [a,ecc,0,inc,o,nu,E,M,dcm,rmag] = param(rl,vl) 

% Constants 
rtd=180/pi; 
dtr=l/rtd; 
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% Compute the argument of perigee 

% Compute the rotation matrix 
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dcm=rot(0,inc,o); 

E = acos( (ecc + cos(nu*dtr) ) / (1 + ecc*cos(nu*dtr)) ); 
M = E - ecc*sin(E) 
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